A systematic approach for developing mechanistic models for realistic simulation of cancer cell motion and deformation

Understanding and predicting metastatic progression and developing novel diagnostic methods can highly benefit from accurate models of the deformability of cancer cells. Spring-based network models of cells can provide a versatile way of integrating deforming cancer cells with other physical and biochemical phenomena, but these models have parameters that need to be accurately identified. In this study we established a systematic method for identifying parameters of spring-network models of cancer cells. We developed a genetic algorithm and coupled it to the fluid–solid interaction model of the cell, immersed in blood plasma or other fluids, to minimize the difference between numerical and experimental data of cell motion and deformation. We used the method to create a validated model for the human lung cancer cell line (H1975), employing existing experimental data of its deformation in a narrow microchannel constriction considering cell-wall friction. Furthermore, using this validated model with accurately identified parameters, we studied the details of motion and deformation of the cancer cell in the microchannel constriction and the effects of flow rates on them. We found that ignoring the viscosity of the cell membrane and the friction between the cell and wall can introduce remarkable errors.


Materials and methods
Lattice Boltzmann method and dissipative coupling. All 32 using its Object-in-Fluid (OIF) module 33 for modelling viscoelastic cells interacting with fluid flow and microchannel walls. Customarily, at the macroscopic level, fluid flow is treated as a continuum and analyzed by solving the Navier-Stokes equations. However, in recent years, LBM became a powerful method of choice, especially at the cellular length scale, due to its high calculation efficiency and its intrinsic parallelism 15 . In LBM, the fluid flow is discretized in time and space by considering the fluid as fictitious particles. In the present work, three-dimensional 19-velocity cube lattice scheme (D3Q19) has been utilized in the LBM governing equations as follows: where n i (x, t) , δ t , δ x , τ and f i (x, t) are the density distribution function, time step, grid spacing, relaxation time toward the equilibrium distribution n i eq , and external force, respectively. At each lattice site, the macroscopic fluid density ρ and velocity u can be obtained from the particle density functions as follows: where, e i is the discrete velocity vector 34 .
Immersed boundary nodes were defined using vertices of triangular mesh representing the surface of the cell membrane. Motions of these nodes due to fluid forces and viscoelastic properties of the cell membrane are governed by the Newton's second law of motion: (1) n i (x + δ x , t + δ t ) = n i (x, t) − 1 τ n i (x, t) − n i eq (x, t) + f i (x, t) for i = 1, 2, . . . , 19 (2) ρ(x, t) = i n i (x, t) for i = 1, 2, . . . , 19 www.nature.com/scientificreports/ where m j is the mass of the node j, and F j is the force applied to that node. The viscoelastic behaviour of the cell is defined by elastic and viscous parameters which will be explained in "Spring-network model of cell". In the dissipative coupling method, the coupling are done by means of the drag force exerted by the fluid on the nodes of the cell membrane 33

as follows
where v is the velocity of the cell-membrane node, and u is the velocity of the fluid at the same position. In this equation, ξ is the phenomenological friction coefficient between fluid and the immersed object that was already calibrated using experimental data 35 . It was reported that ξ is dependent on the surface area of the immersed object ( S ), the number of cell membrane mesh nodes ( n ) and the dynamic viscosity of the fluid ( υ ). Equation 6 defines ξ for a sphere with radius r and discretized into n nodes as follows: where ξ ref is the calibrated friction coefficient for a reference sphere, given in Table 1 Spring-network model of cell. We used the viscoelastic cell model implemented in OIF in which the cell is discretized employing a two-dimensional triangular network of parallel springs and dashpots according to the Kelvin-Voigt viscoelastic model. In this model, the deformation behaviour of the cell is imitated using one viscous damping coefficient and five elastic moduli (stretching, bending, local area, global area, and volume) as described below 15 .
Stretching modulus. The stretch of a spring spanning as an edge on the cell membrane mesh connecting cell membrane vertices A and B is defined as 23 : where k s , L 0 , L , and − → p AB are the stretching modulus, the relaxed length of the edge AB, the current length of the edge AB, and the unit vector pointing from vertex A to vertex B, respectively. Moreover, = L L 0 is the stretch and κ( ) describes the neo-Hookean hyperelastic behaviour of the edge 36 : Bending modulus. Bending is defined in terms of the changes in the angles between two neighbouring triangles. For triangle ABC, the bending force can be calculated as follows 23 : where θ 0 and θ are the angle between two triangles in the resting and current shapes, respectively. k b is the bending modulus and − → n ABC is the unit vector normal to the triangle.
Local area, global area and global volume moduli. The biology of the cell membrane necessitates resisting to changes in the local area, global area and volume of the whole cell 15 . For node A, forces applied on each node of the triangle ABC for satisfying the local area, global area, and global volume preserving constraints can be computed from: www.nature.com/scientificreports/ where S 0 ABC , S ABC are the areas of the triangle ABC in the resting and current shapes, t a , t b , t c are the distances of the triangle vertices to its centroid T. S 0 and S are the areas of the entire cell in resting and current shapes, V 0 and V are the whole volumes of cell in the resting and current shapes. Also, k al , k ag , and k V are the local-area modulus, global-area modulus, and global-volume modulus, respectively. As stated before, − → n ABC is the unit vector normal to the triangle ABC 23 .
Viscous damping coefficient. Since the cell membrane has a lipid layer within its bilayer, it exhibits viscoelastic behaviour. In OIF the membrane viscosity was implemented using a Newtonian viscous damper 15 : where k visc is the viscous damping coefficient, and − → v A is the velocity of node A.
Cell-wall interaction. At a constant pressure drop between the two ends of the constricted microchannel, the passage time (entry time + transit time) of a cell depends on cell deformability as well as cell-wall surface friction. Furthermore, in ESPResSo repulsive forces were defined between the nodes of the cell and microchannel wall to avoid cell penetration into the microchannel wall in order to model the behaviour of the cell near the walls as 15 : where d is the distance between the nodes of cell and wall, d cut is the threshold of repulsive force activation, a and n are the repulsive force coefficients, and − → m is the unit vector pointing from the wall node to the cell node 15 . The repulsive-force parameters were considered to be constant in all simulations in this study.
On the other hand, since originally ESPResSo and OIF did not have friction-force implementation in the source code, to realistically model cell passing through a constriction, we defined the below friction force where µ f is the surface friction coefficient.
Numerical model of cell passage through a microfluidic constriction. In this study, our aim is to develop a validated spring-network model for cancer cells as well as a numerical tool that can accurately find the parameters of this model. The validated model with accurate parameters then can be used to model the deformation behaviour of cancer cells in various scenarios. We identified parameter values that enabled the model to replicate the experimental data 9 of a cancer cell passing through a microfluidic constriction. In their study, a suspended microchannel resonator with a constriction was used to investigate the deformability and surface friction of various types of cancer cells including the human lung cancer cell line H1975. The constriction was 6 μm in width, 15 μm in depth, and 50 μm in length with a 45° tapered entrance. Cells were kept in cell solution (RPMI 1640) at 37 °C and sent through the microchannel with a constant pressure drop that drove the flow in the fluidic microchannel. The flow rate of the blank media at the drop pressure of 1.5 psi was reported to be 38 μL/h. The measurements of the passage of H1975 were performed at two different drop pressures of 1.8 psi and 0.9 psi 9 . The flow rates at these two drop pressures were calculated using Poiseuille's law as listed in Table 2. In this table, the entry time is the required time for the cell to completely enter the constriction and the transit time is the time the cell takes to travel in the constriction and to completely exit it. Because their experimental results 9 were reported based on the cell buoyant mass, to find the cell radius in their experimental data, we used the following equation that relates the cell radius r to its buoyant mass: www.nature.com/scientificreports/ where, m b is the cell buoyant mass, ρ is the cell density, and ρ f is fluid density 9 . We modelled the cancer cell as a sphere with the radius of 6.5 μm and discretized with a triangular surface mesh with 393 nodes. We performed Python scripting in ESPResSo to create the geometrical setup of the microchannel and the computational lattice for our simulations. Table 3 shows the parameters used in our numerical simulations of cancer cell motion and deformation in the microchannel. Because, in the experiments, walls of the microfluidic channels were not lined with endothelial cells, receptorligand interactions between the CTCs and walls did not exist and therefore they were not considered in the present model.
To achieve the goal of this study, a large number of numerical simulations should be executed repeatedly to find accurate model parameters. Since numerical simulation of the whole microchannel setup was computationally expensive, performing the calculations with a smaller model with an acceptable computational accuracy could greatly speed up the parameter identification process. We performed numerical simulations of the motion and deformation of lung cancer cells in four geometrical domains shown in Fig. 1. In this figure, Model #1 (panel a) includes the entire fluid passage, while Models #2, #3 and #4 (panels b, c and d) only include a portion of Model #1. Figure 1e shows the fluid velocity profiles at the depth of 7 μm at the cross-section E-E (shown with a dashed line in Fig. 1b) for the flow rate of 45.6 µL/h. All three partial models, that have smaller geometrical domain compared to the complete microchannel model (Model#1), have velocity profiles very close to that of Model#1. The numerical results of lung cancer cell entry time along with the required simulation runtime were compared in the four models with similar viscoelastic parameters for four various combinations of the parameter values shown in Table 4. The calculated cell entry time and simulation runtime of Models #2, 3 and 4 against those of Model #1 were compared in Table 4. For Models #2, 3 and 4, the average errors were 4%, 4.2% and 1419.2%, respectively, while the average runtimes were 1509.25 s, 267.25 s, and 2161.5 s, respectively. Because of its substantial error, Model #4 cannot be reliably used instead of Model#1. Although Model #2 was 0.2% more accurate than Model #3, its runtime was over five times of that of Model #3. Since Model #3 can sufficiently and efficiently represent Model #1, we used Model #3 in the calculations of the present study. Figure 1f illustrates the fluid velocity profile at the E-E cross section (Fig. 1b) for the flow rate of 22.8 µL/h with different lattice grid spacing for the fluid domain. As Fig. 1f shows the fluid flow is not sensitive to the grid spacing of equal or less than 1.5 μm. Therefore, in this study the grid spacing of 1 μm was used for the fluid domain to have a high accuracy while keeping the computational cost manageable.
The deformation of the spring-network cell model, however, is sensitive to the mesh size so that all parameters should be changed when the mesh size is changed. Therefore, for each mesh size of the cell, the proposed parameter identification process should be repeated.
Genetic algorithm. The spring network model for the cell includes several parameters which need to be set accurately and because the parameters are dependent on one another, finding the best set of parameters to imitate large deformations of cells is a complicated task. In this work, we propose an optimization-based method to accurately identify the parameters of spring-network models of cells. We defined an error function and developed an optimization algorithm based on GA to accurately and automatically identify parameters to minimize the error of the model results compared to the experimental results.
The steps of the GA code are as follows: 1. Random n = 128 bits fractional binary numbers were generated for each model parameter to constitute the initial population for each iteration. The fractional binary numbers were converted to fractional decimal numbers using the following equations: (17) (0.a 0 a 1 . . . a n−1 ) 2 = (a 0 × 2 −1 + a 1 × 2 −2 + · · · + a n−1 × 2 −n ) 10 = q 10 www.nature.com/scientificreports/ where K L , K U are the lower bound and the upper bound given in Table 5. As an example, the following binary fraction for K s with K L = 0.0001 and K U = 10 and n = 128, 0.01111111110110101010100010001101111001 0011110010001110110111101101010110101010100001000011010101101101101010010011110000001101 10 is the decimal K s of 4.994. The initial population includes 40 different chromosomes for each iteration. 2. Parents were selected from the initial population and then crossover were performed for probabilities of higher than 0.8. Mutations were performed for each bit for probabilities of higher than 0.8.    www.nature.com/scientificreports/ 3. Binary numbers were converted to decimal numbers according to the upper and lower bounds of each parameter as provided in Table 5. 4. Numerical simulation of the corresponding experiment was conducted for each set of parameters, and then the outputs were stored and compared with the experimental results. 5. Each iteration had about 400 distinct parameter sets. In order to generate random numbers based on the best results of the previous iteration, 20 best sets of parameters of the previous iteration, were chosen based on the error function. They were stored and added to the initial population of the current iteration. 6. The algorithm stopped if one of the below conditions occurred: (a) The error function value is less than a user-defined threshold.

Validation of the parameter identification method for spring-network model of cells. In order
to assess the sufficiency of our proposed GA-based method for identifying the values of parameters of the spring-network model of cells, we used the method to identify the parameters of the RBC. We used the RBC stretch experiment performed using optical tweezers by Mills et al. 3 . Figure 2a shows the 3-D model of the RBC with 374 mesh nodes that we used in our calculations. The RBC was stretched in the axial direction as shown in Fig. 2b. We used the following error function where s a and s t are the axial and transversal diameter of RBC, respectively, which were calculated numerically at the last time step. Also, e a and e t are the axial and transversal diameter of the RBC measured in stretch experiment with optical tweezers 3 . The stop condition error threshold was defined to be one 15 . Figure 2c compares the experimental results of Mills et al. 3 with the numerical result of our stretching simulation using the parameter set that our proposed method identified (Table 6) by minimizing the error function of Eq. (16). As Fig. 2d shows, the error function greatly decreased from 4.8 to close to 1 after only five iterations and then decreased to below the threshold of one after ten iterations and the calculations stopped and model parameters were identified (Table 6). Figure 2c also shows the numerical result that Cimrak et al. 15 obtained using their manually adjusted parameters given in Table 6 on an RBC with the same number of nodes and mesh configuration as used in our simulation and shown in Fig. 2a. As Fig. 2c and Table 6 show, compared to the parameters obtained by Cimrak et al. 15 , the parameters identified by our GA method enabled closer replication of the experimental results with a substantially lower error (0.73 versus 1.46).
Parameter identification of the spring-network model of cancer cell. In this section, we used our proposed GA method to identify parameters of the human lung cancer cell (H1975). While the time the cell takes to squeeze and completely enter the constriction (entry time) is highly dominated by the deformability characteristics of the cell, the time the cell takes to travel inside the constriction (transit time) is affected by both surface friction and the deformability of the cell 9 . We, therefore, ignored the surface friction effects in calculating the entry time and first identified the viscoelastic parameters during the entry phase by minimizing the following error function: where ET s , and ET e are numerically calculated entry time, and experimentally measured entry time, respectively. Here, n t is equal to 2 for using the entry time at two flow rates.
As described in "Spring-network model of cell", the surface friction force is activated when the distance between the wall and cell nodes is less than d cut . We observed that the number of cell nodes which have the active surface friction is not the same in different flow rates. Hence, the surface friction coefficient was identified for each flow rate separately according to the following error function: www.nature.com/scientificreports/ where i is the experiment number according to Table 2 and TT s , and TT e are numerical transit time and experimental transit time, respectively. Table 7 shows the identified values for the friction coefficients of the lung cancer cell with the radius of 6.5 μm discretized with 393 nodes using experimental results at two different flow rates (given in Table 2). Table 8 shows the entry time, transit time and the calculated errors at these two flow rates.
Steps of cancer cell deformation in constriction. The Table 6 and (3) numerically using the parameters identified using the proposed genetic algorithm identification method. The model with parameters identified using the proposed method can replicate experimental data more closely than the model by Cimrak et al. (d) error minimization evolution of the proposed method. The main drop in the error occurred in the first five iterations and the five next iterations fine tuned the parameters and decreased the error to less than the threshold of one.  Step i), a sudden 22% elongation (Fig. 3a) and a sudden increase of 0.036 in the area strain (Fig. 3c) occurred. The sharp slopes in the leading and trailing edge traces (Fig. 3b) show that the cell moved inside the channel while deforming from a sphere to an ovoid (Fig. 3d).
Between Steps i and viii, the gentle slopes of the leading and trailing edge traces (Fig. 3b) show continuous but sluggish progression of the cell into the constriction. However, details of the deformation are not visible in the leading and trailing edge traces. Pushing the cell against the constriction wall between Steps i and ii (Fig. 3d), caused the cell length to drop to near only 3% longer than the original diameter ( Fig. 3b) but the almost constant area strain shows that the cell expanded in the direction normal to the view of Fig. 3d. Between Steps ii and iii, the cell underwent another rapid elongation (7% increase) and a rapid increase in the area strain (to 0.05) as Fig. 3a,c show. The major portion of the time (89%) during entry was spent between Steps iii and viii when the cell underwent small gradual elongation and gradual increase in the area. In Step vi, the area strain stayed almost constant while the length gradually increased. In Step vii, the cell ceased the constant area regime and eventually completely entered the constriction in Step viii. This was accompanied by a sudden 0.003 decrease in the area strain while the cell suddenly elongated in the narrow constriction to the eventual 29% increase compared to the original diameter.
Step ix is the transit step during which the cell length and area strain remained almost constant while the cell rapidly travelled the narrow constriction as demonstrated by the steep traces of the leading and trailing edges in Fig. 3b.
At the higher flow rate of 45.6 µL/h , the steps described above are completed faster as panels of Fig. 4 show (steps are labeled in Fig. 3). The initial stretch of the cell in Step i was more at the higher flow rate than it was at the low flow rate of 22.8 µL/h as Fig. 4a shows (16.6 vs 15.9 µm). Furthermore, the area strain in Fig. 4c shows a higher jump of 0.07 in Step i than the jump of 0.05 that the cell underwent at the lower flow rate. Between Steps i and viii, unlike at the lower flow rate that the area strain showed the overall continued growth 0.05 to about 0.07, at the high flow rate, the area strain remained almost constant near 0.07 with small undulations (Fig. 4c). Figure 4b shows the axial location of the centre of the cell at both flow rates. As this figure shows, between Steps i and ii, the cell entered further into the narrow constriction at the high flow rate than it did at the low flow rate. This is evident from the cell-centre axial locations of 11.6 µm with the high flow rate versus 8.2 µm with the low flow rate (both at 0.3 ms). This observation is consistent with the cell length plots of Fig. 4a where between Steps i and ii, the cell always elongated more with the high flow rate than it did with the low flow rate. As Fig. 4a shows, while the time span between Steps i and ii at the two flow rates were almost the same, cell deformation evolution between Steps ii and viii were compressed in time at the high flow rate compared to that at the low flow rate. The straight horizontal lines (Step ix) at the end of the plots in Fig. 4a,c show that the transit time was much shorter at the high flow rate than it was at the low flow rate. This is consistent with the almost vertical line at the end of the cell-centre trace of Fig. 4b that indicates very short transit time compared with the steep-sloped end of the cell-centre trace at the low flow rate.

Parsimony of spring-network model for cancer cells.
While the elastic parameters used in the model that we present here for cancer cells are widely used in spring-network models of cells, the viscosity is not included in some cell models 24 . Also, the friction between the cell and the wall is usually ignored in models 9,14 . In the following sections, we investigated how ignoring these two parameters, in order to reduce the number of the parameters of the cancer-cell model, can affect the validity of the model results.

Effects of membrane viscosity on cell deformation in constriction.
Considering viscosity of the cell membrane is crucial in developing a sufficient model of cancer cell motion and deformation. Indeed, the cell deformation in the constrictions showed flow-rate-dependent sensitivity to membrane viscosity as demonstrated in this section. www.nature.com/scientificreports/  (Fig. 5a) and area strain (Fig. 5b) show, although the cell without viscosity reached to Step v (defined in "Steps of cancer cell deformation in constriction") sooner than the cell with the correct viscosity did (5.2 vs 6.2 ms respectively). However, it remained in Step vi (indicated with the straight horizontal line after 5.2 ms) while the cell with the correct viscosity finished all steps and exited the microchannel. Consistently, the 3-D views of the cell in Fig. 6 show that the cell without viscosity remained unchanged at the entrance of the narrow constriction at the shown instances of 6.53 ms and beyond. Meanwhile, the cell with the correct viscosity exhibited a series of deformations that eventually resulted in the entrance to the narrow constriction at 16.53 ms.
On the other hand, at the high flow rate (45.6 µL/h), the cell with ignored viscosity reached Step ii sooner than the cell with the correct viscosity did (1.6 vs 2 ms) as it is clear in the cell length plot of Fig. 5c and the area strain plot of Fig. 5d. The cell with ignored viscosity completed all subsequent steps sooner than the cell with the correct viscosity did and exited the microchannel with the passage time of 1.9 ms compared to that of 2.2 ms for the cell with correct viscosity.
motion of the cell in the microchannel represented by the trajectories of its centre of mass, and (c) evolution of the area strain of the cell (change in the cell area relative to its undeformed shape divided by its undeformed initial area) as the cell advanced in the microchannel. The evolution of the cell at the low flow rate and the steps are marked in Fig. 3 and here it is shown for reference. The cell underwent a higher initial elongation and area strain in Step i with the high flow rate than it did with the low flow rate. The initial high area-strain value, that the cell reached to in Step i with the high flow rate, stayed almost constant during the entry steps. However, with the low flow rate, the area strain gradually increased during the entry steps until it eventually reached to the high value that the cell reached to Step i with the high flow rate. www.nature.com/scientificreports/ Effects of friction on cell passage in constriction. In this section, we demonstrate that in order to develop accurate spring-network models of cancer cells in interaction with a wall, using an accurate friction coefficient is crucial and this parameter should not be ignored. Figure 7 shows the transit step (Step ix as defined in Fig. 3) for the flow rates of 22.8 and 45.6 µL/h at four axial locations along the narrow constriction. In both Fig. 7a,b, the top rows show the cell motion states simulated while ignoring friction and the bottom rows show the cell motion states simulated using the friction coefficient identified for the flow rate (Table 7). We present the location-matched states of cells because the cell in the frictionless condition travelled so fast that time-matched comparison is not useful. The cell is coloured with the velocity magnitude. At the low flow rate (22.8 µL/h ), at the 18 µm axial location, the cell with identified friction still had nonuniform velocity with the maximum speed of about 2 cm/s while the frictionless cell already exhibited rigid body motion with the uniform speed of 3 cm/s. At this flow rate, both with-friction and frictionless cancer cells showed rigid-body motion at the 30 and 42 µm axial locations and the frictionless cell travelled twice as fast as the cell with correct friction did (at the constant speed of about 3 vs 1.5 cm/s). At the higher flow rate (45.6 µL/h ), at the 18 µm axial location, cells in both conditions reached to rigid body and constant-speed motion. The frictionless and with-friction cells travelled at 10 and 7 cm/s, respectively. Ignoring the cell-wall friction caused a remarkable error in the transition time (81% and 55% error for 22.8 and 45.6 µL/h , respectively).
As explained in "Parameter identification of the spring-network model of cancer cell", in the proposed parameter identification method, the friction coefficient should be identified separately for each flow rate. When instead of the friction coefficients individually identified for each flow rate, the average of the friction coefficients identified for each flow rate was used, the following was observed. For the higher flow rate simulation became unstable and crashed, since forces exerted on the cell membrane nodes were high at this flow rate. For the lower flow rate, using the average friction coefficient caused the calculated transit time to have %47 error in comparison with the experimental data.

Discussion
Spring-networks have been increasingly employed for modeling the motion and deformation of red blood cells 24,27,37 , leukocytes 38 , platelets 39 and cancer cells 40 . Spring-network models were already used for designing microfluidic chips 41 and for such designs having a validated model with accurate parameters, as our method provides, is crucial. Furthermore spring-network models can be used to obtain fundamental understandings about metastatic spreading in the body and it is crucially important to use accurate parameters in these models. To the best of our knowledge, our proposed method is the first systematic method to find parameters of springnetwork models to develop validated cancer cell models. www.nature.com/scientificreports/ The parameters of the spring-network model of the RBC were determined by linking a continuum-mechanics model and a spring-network model of the cell 15 . However, that approach needs special assumptions which restrict the model from being general since the derivation is only valid for one specific geometry and cell type. Also, the parameters determined employing that approach do not accurately model the cell behaviour in comparison with experimental results and further calibration of the model parameters is required 37 . Our proposed method directly identifies the spring-network model parameters and does not rely on any other intermediate models or other restrictive assumptions.
To assess the effectiveness of the proposed method, we applied it to find parameters of the spring-network model of the RBC to imitate its deformation in the stretch experiments performed using optical tweezers. We showed that the parameters obtained using the proposed method enable closer replication of the experimental data compared to the reported parameters in the literature that were adjusted manually.
Ye et al. 42 proposed functions that associated the transit time with bending modulus alone or shear modulus alone while treating all other parameters fixed. They later expanded the function to include both bending and shear moduli 41 . They studied effects of these two moduli on the entry time and fitted functions to relate the two moduli with transit time. Our approach in this work is different as we proposed a method to identify model parameters based on experimental data of cell motion in a microfluidic constriction. www.nature.com/scientificreports/ The model of the cancer cell used in this study only includes the cell membrane represented by a network of springs. However, the deformation of the cell is also affected by the elasticity of the cell nucleus and cytoplasm. We identified the parameter values of the spring-network model of the cell using experimental data of the deformation of the entire cell. Therefore, the elasticity of the cell nucleus, cytoplasm and membrane are all lumped together and are captured by the parameters of the spring network of the cell.
For the two flow rates for which we had experimental data, the discrete friction coefficient was found to strongly depend on the flow rate. To correctly capture the physics of the interaction of the cell with the channel wall, in terms of the cell transit time, the friction coefficient was found to decrease from 0.0118 at 22.8 µL/h to 0.0009 at 45.6 µL/h. We observed that, compared to the low flow rate, at the high flow rate the cell nodes got closer to the channel wall causing a higher nodal normal repulsive force. The proposed algorithm balanced the intense normal repulsive force with a reduction in the coefficient of friction at the high flow rate (Eq. 15). These observations, however, were based on the experimental data at only two flow rates and should be verified by collecting experimental data at multiple flow rates and performing the parameter identification using the method proposed in this work.
An immediate application of the proposed method with the accurately identified parameters of the springnetwork model of the cancer cell is to predict the cell deformation and passage in the microchannel constriction at flow rates not measured experimentally. In the absence of experimental data to identify the friction coefficient for flow rates other than the two flow rates, a linear relationship between the friction coefficient and flow rate was assumed using µ 1 f and µ 2 f in Table 7. Figure 8a,b illustrate the entry time and transit time of the cell as functions of flow rate, respectively. As Fig. 8a shows, the cell entry time to the narrow constriction generally follows a power law decreasing trend with increasing the flow rate from 15.5 ms at 23 µL/h to 1.6 ms at 46 µL/h. Similarly, the general trend of cancer-cell transit time in the narrow constriction is a power law decrease with increasing the Figure 7. Effects of friction between the walls of the microchannel and the cancer cell on the cell motion and deformation during its transit in the constriction (Step ix in Fig. 3) of motion at the flow rates of 22.8 µL/h (a) and 45.6 µL/h (b). The cell is viewed from the top (perpendicular to the microchannel axis). Because neglecting friction caused the cell to pass very fast, instead of time-matched states, axial location-matched states are shown. At both flow rates, the cell velocity magnitude at any location in the constriction was always overestimated when the friction was ignored.  (Fig. 8b). Interestingly, as Fig. 8c shows, the transit time and entry time have a linear relationship with a coefficient of determination of 0.92. In this figure the increasing directions of both axes correspond to decreasing the flow rate as Fig. 8a,b show. The values of the flow rate, entry time and transit time reported in Fig. 8a,b are important for designing new experiments for obtaining cell-deformation information at several flow rates for further improving the accuracy of the identification of the parameters of the spring-network model of cancer cells. Cell deformability was shown to be a promising label-free biomarker for cancer diagnosis 43 and some high-throughput single-cell microfluidic devices have been proposed for cell characterization based on this biomarker 7,44,45 . Moreover, recently artificial intelligence was used to substantially increase the specificity of deformability-based cell sorting in microfluidic chips 46 . The method proposed in the present study can be used to obtain validated spring-network models of motion and deformation of cells that can be used for simulating passage of cells in microfluidic devices. This enables simulation of the cell deformation for designing new microfluidic devices by changing and optimizing the geometrical and fluid flow parameters before manufacturing the chip for cell characterization based on deformability.
Imaging modalities such as magnetic resonance phase contrast angiography in combination with mathematical optimization can be used for obtaining patient-specific vascular networks 47 . The validated cancer cell models  www.nature.com/scientificreports/ also enable computationally efficient simulations of the metastatic processes in such vascular networks. These models can include cell-wall friction that can affect CTC adhesion and extravasation. A sufficient model of deformation and motion of cancer cells should include viscosity of the cell as we demonstrated that ignoring the viscosity of the lung cancer cell had remarkable effects on its motion and deformation in the constriction. At the low flow rate (22.8 µL/h), the membrane viscosity encouraged deformation so that the cell completely squeezed in the microchannel sooner than it did in the absence of viscosity. However, at the high flow rate (45.6 µL/h), the cell viscosity resisted cell deformation and slowed down the deformation and therefore the entry time with viscosity was longer than it was without viscosity. This shows the sophisticated effects of the cell-membrane viscosity on the deformation and motion of the cancer cell and that these effects cannot be numerically achieved with a simplistic model.
A sufficient model of a cancer cell passing through a constriction should also consider the friction between the cell and the microfluidic wall as we demonstrated in this study. We showed that for a spring-network model, the flow-rate dependent friction coefficient should be identified and used in the model to accurately simulate the interactions of the cancer cell with the microchannel wall.
The method proposed in this study, the validated model and the knowledge obtained from it can be the basis for models of CTC in circulation, entrapment and extravasation stages of metastasis. Such models are important for obtaining better fundamental understandings of the metastasis process as well as for developing novel predictive tools.

Limitations.
Higher accuracy in the simulation may be attained if the model parameters are obtained for the cell with a higher number of nodes. However, this would increase the computational cost to the extent that it cannot be possible for the GA code to find an accurate solution for the model parameters based on the experimentally measured passage time. The algorithm proposed in this work should be improved to enable more computationally efficient parameter identification process for finer discretization.
Furthermore, as it has been shown for parameter identification of the RBC using the stretch experimental data, the GA code can efficiently and accurately identify model parameters when the number of inputs and outputs are increased. For cancer cell deformation in a microchannel, measuring the output (e.g., entry time) at several flow rates can increase the number of input-output pairs and therefore can enable more efficient and more accurate identification of the cancer-cell model parameters compared to what we reported here. Future experiments should be carried out considering this requirement.
In this study only one lung cancer cell line (H1975) was investigated. Future studied should consider other cell types and the performance of the proposed method for those cell types should be investigated.
This study only considered flow rates of up to 45.6 µL/h. Future experimental and modelling studies should consider the effects of higher flow rates on the modes of motion and deformation of the cell in the constriction and should investigate whether friction and membrane viscosity parameters still play crucial roles at such high flow rates.
For the two flow rates for which we had experimental data, the discrete friction coefficient was found to strongly depend on the flow rate. To correctly capture the physics of the interaction of the cell with the channel wall, in terms of the cell transit time, the friction coefficient was found to decrease from 0.0118 at 22.8 µL/h to 0.0009 at 45.6 µL/h. We observed that, compared to the low flow rate, at the high flow rate the cell nodes got closer to the channel wall causing a higher nodal normal repulsive force. The proposed algorithm balanced the intense normal repulsive force with a reduction in the coefficient of friction at the high flow rate (Eq. 15). These observations, however, were based on the experimental data at only two flow rates and should be verified by collecting experimental data at multiple flow rates and performing the parameter identification using the method proposed in this work.