Research Article  Open Access
Jun Qiu, FangFang Li, "Identification of Viscoelastic Constitutive Parameters of a Cell Based on FluidStructure Coupled Finite Element Model and Experiment", Mathematical Problems in Engineering, vol. 2019, Article ID 5868561, 13 pages, 2019. https://doi.org/10.1155/2019/5868561
Identification of Viscoelastic Constitutive Parameters of a Cell Based on FluidStructure Coupled Finite Element Model and Experiment
Abstract
Although mechanical properties of a single cell as well as its mechanical response to stimulation are significant for mechanotransduction, it is difficult to quantitatively identify the viscoelasticity of a cell. This study proposes a methodology to determine the viscoelastic parameters for a single cell combining flow chamber experiment and fluidstructure coupled finite element model. The observation from the experiment is used to compare with the deformation response of the cell with different parameters in the finite element model with standard linear solid constitute model, which are adjusted until the two displacement response curves from the experiment and the simulation accord with each other, and then the corresponding combination of the viscoelastic parameters is determined to represent the viscoelasticity of the cell. The proposed methodology is applied on osteocytes in this study but can be generalized to other cells. The results indicate that both and of an osteocyte are in the order of several hundred Pa and is in the order of kPa·s.
1. Introduction
Establishing and improving mechanical models of a cell is an active research field in biomechanics [1]. Various cellular mechanics models have been proposed, mainly including two categories: the models based on micro/nanostructure and the continuum cell models [1–3].
The models based on micro/nanostructure consider the cytoskeleton as the main structural component [2–9]. The mechanics of cytoskeleton of adherent cells were studied [3]. For suspended cells, such as red blood cells, a model of microhaemoglobin network [10, 11] was established to study the role of cell membrane and haemoglobin network in erythrocyte macromorphism.
On the other hand, the continuum methods regard a cell as continuums with determinate material properties. The parameters based on constitutive models of certain materials can be measured by experiments. Although these models do not provide an indepth study of the mechanics between biomolecules and bioions, the continuum methods are effective in understanding the overall mechanical properties of cells at the cellular level. Once the continuum mechanical model of a cell is established, the details of stress and strain distribution on the cell can be provided, which helps to understand the distribution and transmission of force between cytoskeleton and subcellular components. The cell mechanics models based on the continuum assumption [12] include solid model, mainly used for adherent cell research [13]; cortical coating droplet model, mainly used for suspension cell research [14]; structural damping model with powerlaw, mainly used for the study of dynamic properties of adherent cells [15]; biphasic model, mainly used for the study of muscle and skeletal cells [16]. Other cell mechanics models are refined based on these basic models.
In experiments, the cell size is between ten microns and dozens of microns, and the thickness of cell membrane is between several nanometers and dozens of nanometers [13]. Conventional macromechanical loading methods and experimental techniques cannot be directly used in the study of cell mechanics. With the advancement of micromanipulation technology, many mechanical loading experimental methods for cultured cells in vitro have been proposed, and different types and functions of experimental devices have been developed to achieve quantitative study of the mechanical properties of cells, including micropipette aspiration [14, 17], unconfined cell compression [18], cytoindentation [19], atomic force microscopy (AFM) [20, 21], magnetic bead rheometry [22], and optical traps [23]. In particular, many researchers have studied the adhesion and deformation of adherent cells under shear flow [24], for example, the shear effect of blood on endothelial cells in vivo [7] and the adhesion and debonding of cultured cells under shear stress [24]. The stress and deformation of a cell have direct impacts on the function and structure of the cell. Because of the extremely complex cytoskeleton network system inside the cell, the cell has the ability to deform actively and resist passive deformation. The mechanical properties of a cell are closely related to its structure and function. Viscoelasticity is an important part of in the research of cellular mechanical properties.
Osteocytes, mature bone cells embedded in mineralized bone matrix, are under dynamic shear fluid in vivo and exhibit significant viscoelastic behavior [25]. Abundant evidence has shown that osteocytes are the key mechanosensor cells that directly regulate boneforming osteoblast and boneremoving osteoclast activities [26, 27]. Thus, osteocytes are critical to etiology and new treatments for osteoporosis [28].
In this study, the viscoelastic parameters of osteocytes are determined by a comparison between finite element model and experimental observation. A flow chamber experiment is carried out to record the behavior of osteocytes under shear flow, which is simulated by a fluidstructure interactive finite element model. Standard linear solid (SLS) model is used to describe the viscoelastic behavior of cells, the three parameters in which are determined by comparing the deformation of the cell under shear flow stress. Applying the methodology on three osteocytes, the three parameters of viscoelasticity are identified. The results show that both and are in the order of several hundred Pa and is in the order of kPa·s.
2. Methodology
2.1. Flow Chamber Experiment
In cell flow chamber experiment, osteocytes adhere to the basal glass plate, which is placed in a square flow chamber, and a unidirectional flow of cell culture medium flows through the chamber. Osteocytes are deformed by fluid shear force, as shown in Figure 1(a).
(a)
(b)
(c)
The length of the flow chamber is 5 cm, which is three orders of magnitude larger than the cell scale. The height of the flow chamber is 550 μm and the width is 700 μm. The size of osteocytes is around 10 μm.
It is generally believed that the shear stress of adherent cells in microvessels is 10dyn/. Based on the theoretical analytical relationship between wall shear stress and flow rate, the flow rate in a square cavity can be estimated. In the experiment in this study, the flow rate of 1.318 ml/min is strictly controlled by using a microinjection pump, so as to control the shear force of flow field on cells.
In the cell shear flow experiment, the bottom and side deformations of the cells are obtained by a selfbuilt quasithreedimensional microscopic imaging system, as illustrated in [29]. The microscopic images of the cell are recorded at 12 Hz frequency, and the displacement information of any point on the cell can be obtained by digital image correlation algorithm [29]. Figures 1(b) and 1(c) show the side images of the osteocyte before and after the deformation under shear stress.
2.2. D FluidStructure Coupled Finite Element Model
Software ADINA developed by the research group led by K. J. Bathe in 1975 is a dynamic nonlinear finite element analysis software. ADINA has remarkable ability to solve fluidstructure interaction problems. In the process of fluidsolid coupling calculation, the fluid force is applied to the structure through the interface of fluid and solid, and the deformation of the structure affects the boundary conditions of the fluid in turn. When doing fluidstructure coupling calculation, the structure model and fluid model need to be built in ADINA structure module (ADINA Structures) and fluid module (ADINA CFD), respectively, and then the derived finite element input files from the two models are put into the ADINA fluidstructure coupling solver (ADINAFSI) for calculation. Complex multifield coupled physical problems can be simulated and predicted by ADINA coupling solution, which is adopted in this study.
In order to verify the computational accuracy of the fluidsolid coupling calculation program and to provide an optimal grid size reference for 3D modeling, we studied the 2D model before establishing the 3D model.
2.2.1. Basic Assumptions of Finite Element Model
① Osteocytes are regarded as homogeneous, isotropic, and nearly incompressible viscoelastic continuum.
② Active force produced by osteocytes under fluid shear is ignored; i.e., only the passive deformation of osteocytes is analyzed.
③ The hydrodynamic properties of water at room temperature are equivalent to those of the cell culture medium. Since the cell shear flow experiment is carried out in the laboratory and the experiment is executed within a few minutes, the influence of environmental temperature on the experiment is not considered.
2.2.2. D Finite Element Model of Cell Shear Flow Interaction
The problem of flow in a square crosssection tube is the basis of the interaction of cell shear flow. The numerical simulation of the flow in a square crosssection tube needs to be realized firstly to study the interaction of cell shear flow. Within the laminar flow range, the fluid is assumed to be incompressible. The crosssection height of the square tube is H and the width is W. For steady laminar flow of incompressible fluid, the NavierStokes equation under unidirectional flow in a square crosssection tube is shown in [24]where is the fluid viscous coefficient, is the fluid pressure, is the flow velocity, is the dimension of the flow field direction, is the direction of the height of the flow chamber, and is the direction of the width of the flow chamber. The boundary conditions for this problem are shown in where is the height of the flow chamber and is the width of the flow chamber.
When the pressure gradient at the flow direction, the geometrical size of the square tube, and the viscous coefficient of the fluid in the tube are known, the distribution function of the velocity profile under the laminar flow state can be obtained as in
The wall shear stress is as in
The flow rate in the square section is as in
In the 2D finite element model of tube flow, the tube thickness is 550 μm and the length is 1400 μm. The wall nonslip boundary conditions are applied along the flow direction. The entrance boundary conditions are pressure P=1 , and the exit boundary conditions are free. Three mesh sizes are designed, including coarse meshes of 18.4 μm, general meshes of 9.2 μm, and fine meshes of 4.6 μm.
Table 1 shows the convergence analysis results of 2D fluid meshes density. It can be seen that the solution of the 2D finite element model of tube flow converges to the theoretical solution by refining the mesh density. Even when the mesh density is small, the flow rate and velocity have high accuracy. Since the shear force is related to the first derivative of the flow rate, it is difficult to achieve high accuracy. Only refining the meshes can derive a more accurate shear force.

Comparing the analytical solutions and numerical solutions, it comes to a conclusion that the results of finite element analysis for such tube flow problems using ADINA software are reliable.
2.2.3. Convergence Analysis of 2D FluidStructure Coupled Finite Element Model
The fluidstructure coupling finite element model of 2D cell shear flow interaction is shown in Figure 2. In the flow chamber, the fluid in contact with cells is discretized by a fine mesh to keep the mesh size of the flow field in accordance with the cell mesh size. In this model, a cell is represented by a semicircle, and the bottom of the cell is completely fixed on the base of the flow chamber; i.e., cell debonding is out of consideration.
The cell and its adjacent fluids adopt three different mesh sizes: 0.4μm, 0.8μm, and 1.6 μm. The responses of different mechanical variables with different mesh sizes are extracted, including the displacement of the cell's highest point along the flow field direction (y direction), the displacement of the cell's highest point along the flow field vertical direction (z direction), the positive pressure distribution along the cell surface, and the fluid shear stress distribution along the cell surface.
2.3. D FluidStructure Coupled Finite Element Model
2.3.1. Fluid Model
In the finite element simulation, the establishment of such fullscale model of the flow chamber requires enormous computational resources without necessity. The length of the flow chamber in the finite element model needs to be determined by theoretical deduction to ensure that, under certain boundary conditions, the fluid can achieve fully developed laminar flow in the chamber.
The Reynolds number is as in
In the cell flow chamber experiment, the flow velocity in the tube can be estimated by the crosssection area and flow rate of the square tube, and the width of the square tube can be taken as the characteristic size , while the density of water is constant. It can be estimated that the Reynolds number of this flow problem is 30. Generally, laminar flow is recognized when the Reynolds number is less than 2100. For laminar flow, there is a prediction formula to estimate the entrance length required for the full development of the flow field, as shown in
In this study, d is 550 μm, and thus the length of the entrance needed for the full development of the flow field is 1400 μm.
According to hydrodynamic experience, if the profile shape of velocity boundary condition applied at the entrance approximates the fully developed velocity profile, the velocity profile at the entrance can evolve to a fully developed velocity profile after a very short distance. It is of great significance to reduce the fluid field in the finite element model and thus to save computational resources and improve computational efficiency. The shape and size of the fully developed section can be obtained from the results of the finite element model of the tube flow with a length of 1400 μm. Applying such fully developed velocity profile as a boundary condition to a shorter (200 μm) flow chamber model, a reliable fully developed flow field near the cell with a smaller computational model can be obtained.
The wall condition is applied around the square flow chamber; the pressure boundary condition (1400 μm flow chamber model) and the parabolic velocity boundary condition (200 μm flow chamber model) are applied at the entrance; and the free boundary is set at the exit.
When determining the fluid domain, the space occupied by cells needs to be removed. A fluidstructure coupling boundary condition is applied to the interface between fluid and cell.
2.3.2. Structure Model
Before establishing the final fluidstructure coupling finite element model based on 3D reconstructed cell configuration, a model based on hemispherical cell configuration is established to verify the validity of the model with simple geometric configuration.
Based on the fluidstructure coupling model of hemispherical cell configuration, the fluid domain is shown in Figure 3(a). The length of the square crosssection tube is 1400 μm. Wall shear conditions are applied around the tube, uniform velocity boundary conditions are applied at the entrance (57.06 mm/s), and free boundary conditions are applied at the exit. It is generally believed that the shear stress suffered by adherent cells in microvessels is 10 dyn/cm^{2}, and then the flow rate in the tube can be estimated according to the theoretical analytical relationship between the wall shear force and the flow rate. The cell is fixed in the middle of the lower reaches of the whole tube. Transition mesh density technique is used to divide the spatial mesh of fluid domain. The mesh density near the cell is high, while the mesh density far from the cell is low. The viscosity of water at room temperature is regarded as the viscosity of the culture medium. The cell mesh and the fluid mesh near the cell are densified, as shown in Figure 3(b). The cell mesh size is 0.4 μm.
On the basis of the model above, the initial configuration of the cell is modified. The experimental data on the real osteocyte is applied to reconstruct 3D cell configuration, as illustrated in [30]. To improve computational efficiency and ensure computational accuracy, the length of the flow chamber is reduced to 200 μm. At the same time, the boundary condition of the velocity at flow chamber entrance is set as a paraboloid, the parameters of which are determined by the calculation results of the extracellular flow field velocity. In the fluid field, the area near the cell needs to be meshed more intensively. The area far from the cell can be divided into larger meshes. The meshes in the fluid domain near the cell can be divided into transitional meshes.
2.4. Viscoelastic Constitutive Model
Previous studies have shown that adherent cells present certain viscoelastic properties [22, 31–33]. The standard linear solid (SLS) can be used to describe the viscoelastic behavior of cells [34, 35], as shown in Figure 3(c).
There are three parameters in a SLS. Establishing the stressstrain relationship in time domain and applying the boundary conditions under stress relaxation conditions, the stressstrain relationship of viscoelastic materials in Laplace space can be obtained by Laplace transformation, the simple expression of which is shown in where is the relaxation modulus of the viscoelastic materials, is creep stress, and is initial strain. It can be seen from the expression of relaxation modulus that, at the initial moment t=0, the elastic modulus is , which is called the instantaneous modulus, while when the time lasts for a long time, the second term of (8) becomes 0; i.e., the elastic modulus becomes , which is called the equilibrium modulus. The viscous coefficient and determine the relaxation response characteristic time of the viscoelastic materials together, which describes the evolution speed of viscoelastic materials from initial deformation to final equilibrium.
The viscoelastic material model in ADINA software is expressed by general viscoelastic constitutive relationship, as shown in
Equation (9) and (10) show the relationship between shear stress and shear modulus, and the relationship between normal stress and bulk modulus, respectively. The shear modulus and bulk modulus are expressed in Prony series, respectively, as shown in
The parameters in Prony series characterize the mechanical properties of viscoelastic materials. Assuming that cells are nearly incompressible materials, their bulk modulus can be regarded to be very large. Thus only the shear modulus of cells needs to be determined in the model, which can be magnified by several orders of magnitude to represent the bulk modulus of cells.
Retaining the constant term and the first exponential term of the shear stress in the form of Prony series, the first two forms of Prony series are the same as the viscoelastic material deduced before. Through the simple equivalence between parameters, the relationship between parameters in Prony series and those in SLS model can be easily found, as shown in
To verify (13), the relaxation problem of viscoelastic materials under uniaxial tension is studied by theoretical and numerical methods, respectively. The theoretical solution is shown in
According to (13), a set of SLS parameters are transformed into parameters in Prony series and input into ADINA model. And for the same relaxation problem, the theoretical and numerical solutions are obtained, respectively, which are in good agreement as shown in Figure 4. Thus the relationship in (13) is valid.
2.5. Identification of Viscoelastic Parameters of a Cell
A fluidstructure coupled finite element model based on 3D cell configuration is established. Applying the same boundary conditions to the model as in the experiment and assigning certain values to the parameters of the cell viscoelastic model , the mechanical response of an osteocyte under the selected material parameters can be calculated. If the mechanical response is the same as that of the cells measured in the experiment, the selected parameters can be considered to be able to characterize the mechanical properties of the osteocyte.
The overall framework of the method determining the viscoelastic parameters is shown in Figure 5. After establishing the fluidstructure coupling model for cell shear flow interaction, a set of estimated values of the three parameters osteocyte is given directly and then input into the finite element model for calculation. These parameters determine the 3D configuration of the cell under shear stress in the fluidstructure coupled finite element model. Meanwhile, the deformation process of cells under the fluid is observed by a selfbuilt quasi3D cell microscopic system, as illustrated in [29]. Quantitative data of cell deformation can be obtained by directly processing the cell images obtained in the experiment. Comparing the cell deformation configurations obtained from the finite element model based on the estimated values of three parameters with the experimental data, if they agree with each other well, the estimated values are taken as the viscoelastic parameters of the osteocyte, while if there is still a large gap between them, the viscoelastic coefficients will be reestimated until the proper values of the three parameters are determine.
Through literature research, we have known the general range of mechanical parameters of osteocytes, which is below 1kPa for suspended or partially adhered osteocytes [36, 37]. That is to say, we choose the initial value by prior knowledge.
After applying the methodology on three osteocytes, three viscoelastic parameters of each cell were obtained.
2.6. Calibration of FluidStructure Coupled Finite Element Model
2.6.1. Calibration of 2D Model
Table 2 shows the calculated values of the fluidstructure interaction model under three groups of meshing. It can be seen that shear stress and pressure depend more on the meshes than cell displacement. Such phenomenon is easy to understand in numerical calculation. Since stress is related to the first derivative of the flow velocity, the stress loses the firstorder accuracy in the process of deriving the velocity. Therefore, it is necessary to study the influence of finer meshes on the shear stress and normal stress for the model.
 
Note: refers to the relative error between the calculated value under the mesh density and that under the adjacent fine meshes; and all the physical quantities are from the highest point of the cell. 
As can be seen from Figure 6(a), the pressure distribution on the cell surface remains basically constant as the mesh size decreases to less than 0.4 μm; i.e., the pressure response on the cell surface has reached a convergent solution.
(a)
(b)
As shown in Figure 6(b), the shear stress in the central part of the cell increases slightly and that in the bilateral part of the cell decreases slightly as the mesh size decreases from 1.6 μm to 0.1 μm. When the cell size is less than 0.4 μm, the distribution of shear stress does not change except for the small area at the top of the cell, which means that the shear stress on the cell surface basically converges. The distribution of shear stress at the top of cells has always been dependent on the scale of the mesh, which may be related to the inflection point of the wall boundary conditions. It is believed that it is no longer meaningful to reduce the mesh to capture the convergent solution of the shear stress. The stress and the corresponding deformation of the whole cell are concerned primarily, so it can be regarded that the shear stress distribution on the cell surface converges when the cell size is less than 0.4 μm.
In Table 3, the stability of the fluidstructure interaction model in calculating the pressure and shear stress distribution along the cell surface with different mesh sizes (0.1μm, 0.2μm, 0.4μm, 0.8μm, and 1.6μm) is summarized. The calculation shows that the shear stress is relatively dependent more on the mesh, and its convergence rate is relatively slow with the decrease of the mesh size. In view of the whole cell deformation, rather than the local stress on cell surface, it can be concluded that when the mesh size is 0.4 μm, the pressure and shear stress on cell surface converge.
 
Note: refers to the relative error between the calculated value under the mesh density and that under the adjacent fine meshes; and all the physical quantities are from the highest point of the cell. 
2.6.2. Calibration of 3D Model
Inputting a set of three viscoelastic parameters of cell viscoelasticity into the fluidstructure coupling finite element model of cell shear flow interaction, the cell deformation information can be obtained, as shown in Figures 7(b)–7(g), in which colors represents the displacement information of cells along the direction of flow. The redder the color is, the larger the displacement is. It can be seen that the displacement at the top of the cell along the flow direction is the most significant, which is consistent with the experimental results.
The top area of the cell is chosen to represent the deformation of the whole cell configuration since it is the most obvious area of deformation and the displacement information here fully reflects the characteristics of the whole cell deformation. Figures 7(a) and 7(b) show the microscopic image of this area obtained in the experiment and the calculation results from the finite element model. Besides, the displacement of the cells along the flow direction is significant under the shear stress, while the displacement along the vertical direction of the flow field is small, so only the displacement of the cells along the flow direction is taken into account.
The relationship between the relaxation modulus of viscoelastic materials along the time is shown in (8).
In the following study of the three viscoelastic parameters of the cell, the range of the parameters found in previous studies are taken as [29, 38]. When adjusting , the other two parameters are kept constant, = 600 Pa and η=1000Pa·s. The displacement at the top of the cell along the flow direction is shown in Figure 8(a). In the second half of the observed response curve, it can be found that significantly affects the response of the cell after equilibrium state, and slightly affects the initial deformation response at the starting point of the curve. Thus the equilibrium response of the cell can be reduces by increasing .
(a)
(b)
(c)
Adjusting the value of and keeping the other two parameters constant, =120 Pa and =1000 Pa·s, the displacement at the top of the cell along the flow direction is extracted, as shown in Figure 8(b). At the starting point of the response curve, it can be found that significantly affects the initial deformation response of the cell. From the second half of the response curve, it can be seen that basically does not affect the response of the cell after equilibrium state. Thus the initial deformation response can be reduced by increasing .
Adjusting the value of the viscosity coefficient and keeping the other two parameters constant, =120 Pa and =600 Pa, the displacement at the top of the cell along the flow direction is extracted, as shown in Figure 8(c). At the starting point of the response curve, it can be found that does not affect the initial deformation response of the cell. Then it can be observed that, with the increase of the viscosity coefficient , the velocity of the cell which tends to equilibrium state slows down and the time needed to reach equilibrium is prolonged.
Through the study above, the effects of each parameter of the three viscoelastic parameters on the cell mechanical response under shear flow are understood, which helps to estimate the increment of each parameter according to the mechanical response of the cell and then to determine the optimal combination of the three viscoelastic parameters efficiently.
3. Results and Discussion
A specific osteocyte is selected to extract the displacement at the top area of the cell along the flow direction from the experimental results, as shown in Figure 9(a). The displacement response is typical as viscoelastic materials, which is also similar to the deformation response curve observed by other researchers [38].
(a)
(b)
(c)
(d)
According to previous reports on cell viscoelasticity, the first group of viscoelastic parameters called M1 is estimated firstly. The mechanical response of the cell represented by M1 in the fluidstructure coupled finite element model is shown in Figure 9(b). It can be seen that there is a certain gap between the mechanical response of the cell based on M1 and that of the cell in the experiment.
When estimating the viscoelastic parameters of the cell for the second and the third time, the values of and are reduced to a greater extent, and the value of is increased slightly. Finally, the mechanical response of the cells represented by the selected viscoelastic parameter combination M3 is in good agreement with that of the cells in the experiment, as shown in Figure 9(c). Therefore, it is considered that the parameters of M3 represent the viscoelastic properties of the cell well, which are the viscoelastic parameters for the osteocyte numbered #01.
The process above shows the determination method for viscoelastic parameters of a particular cell. The actual process has undergone more iterations, and only three representative steps are presented here.
Even for the same kind of osteocyte, each osteocyte has different mechanical properties and responds differently to the same external mechanical stimuli. But interestingly, the adherent osteocytes show viscoelastic characteristics of delayed response to load. Figure 9(d) shows the response of three different osteocytes to the same fluid shear stress in the experiment, where the vertical coordinates of the curve indicate the displacement at the top area of the cell along the flow direction.
Three osteocytes are analyzed, and the final simulative and experimental results are shown in Figure 9(d). The identification results of the three parameters of viscoelasticity corresponding to these three osteocytes are shown in Table 4. Among the three parameters of osteocyte viscoelasticity, both and are in the order of several hundred Pa and η is in the order of kPa·s.

Our results are consistent with previous understanding of osteocytes elastic parameters. Many scholars believe that the stiffness of suspended or partially adhered osteocytes is below 1kPa, while that of fully adhered osteocytes is above 1kPa [36, 37]. The difference between our study and traditional studies on the mechanical properties of osteocytes is that we not only measured the elastic parameters, but also obtained the viscosity parameters of the cells. We believe that the viscosity parameters of cells might be an important indicator of cell mechanics and mechanotransduction.
Some scholars may wonder why a simpler linear viscoelastic constitutive model is chosen to study the viscoelastic behavior of cells here. When we carried out the research, we also tried different viscoelastic constitutive models to describe the cell mechanics properties. After going through lots of literatures, we decided to adopt a linear viscoelastic model with only three mechanical parameters, because the form of linear viscoelastic model is simple, the mechanical significance of each parameter is clear, and the three parameters can be measured through well design experiments one by one. Although more complex nonlinear viscoelastic models may describe viscoelastic behavior more realistically, the corresponding meaning of each parameter is not easy to recognize, nor is it easy to determine the parameters through experiments. Therefore, we finally choose the linear viscoelastic constitutive model in this study.
Instead of using the current popular optimization algorithms, such as genetic algorithm, ant colony algorithm, and other global optimal algorithms, we adopted the basic idea of optimization to find the viscoelastic mechanical parameters of cells. The reasons come from three aspects. First, we have mastered the mechanical significance of the three parameters in the linear viscoelastic model and can predict the changes in the mechanical behavior of cells caused by the changes in parameters. Secondly, through literature research, we have known the general range of mechanical parameters of osteocytes. We only need to search within this range, and the optimal solution found could be considered as the global optimal solution. Thirdly, the 3D fluidsolid coupling model requires a large amount of calculation, and optimization methods such as genetic algorithm and ant colony algorithm need tens of thousands of iterations to get better optimization results, so the calculation cost is too high. Based on the above three reasons, we choose the optimization algorithm based on manual judgment of optimization direction and step size in this paper.
4. Conclusion
A methodology to identify the viscoelasticity of a cell is proposed in this study. A flow chamber experiment is firstly carried out, and the deformation of the cell under shear flow is observed by a quasi3D microscopy system. Meanwhile, a fluidstructure coupled finite element model is established using SLS assumption. The deformation of the cell simulated by the FE model is compared with the observation in the experiment. The three viscoelastic parameters are adjusted based on the comparison between the displace response of the cell in the fluidstructure coupled finite element model and that from the experimental results until the simulative and experimental results accord with each other, when the parameters in the FE model can be regarded to be able to represent the viscoelasticity of the cell. The methodology is applied on osteocytes, which have been hypothesized as mechanosensors in bone. The results indicate that both and of an osteocyte are in the order of several hundred Pa and is in the order of kPa·s. The methodology proposed in this study can be generalized to identify the viscoelasticity of any other cells. The methodology for measuring the mechanical properties of cells developed in this paper should be developed towards the aspect of rapid and realtime measurement, which may contribute to the disease detection methods with cell mechanical parameters as indicators in the future.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by National Natural Science Foundation of China (Grants nos. 51879137, 91847302, and 11402136).
References
 C. T. Lim, E. H. Zhou, and S. T. Quek, “Mechanical models for living cells  A review,” Journal of Biomechanics, vol. 39, no. 2, pp. 195–216, 2006. View at: Publisher Site  Google Scholar
 R. L. Satcher Jr. and C. F. Dewey Jr., “Theoretical estimates of mechanical properties of the endothelial cell cytoskeleton,” Biophysical Journal, vol. 71, no. 1, pp. 109–118, 1996. View at: Publisher Site  Google Scholar
 D. Stamenović and M. F. Coughlin, “The role of prestress and architecture of the cytoskeleton and deformability of cytoskeletal filaments in mechanics of adherent cells: A quantitative analysis,” Journal of Theoretical Biology, vol. 201, no. 1, pp. 63–74, 1999. View at: Publisher Site  Google Scholar
 U. G. Hofmann, C. Rotsch, W. J. Parak, and M. RadMacHer, “Investigating the cytoskeleton of chicken cardiocytes with the atomic force microscope,” Journal of Structural Biology, vol. 119, no. 2, pp. 84–91, 1997. View at: Publisher Site  Google Scholar
 P. Cañadas, V. M. Laurent, C. Oddou, D. Isabey, and S. Wendling, “A cellular tensegrity model to analyse the structural viscoelasticity of the cytoskeleton,” Journal of Theoretical Biology, vol. 218, no. 2, pp. 155–173, 2002. View at: Publisher Site  Google Scholar  MathSciNet
 J. H.C. Wang, “Substrate deformation determines actin cytoskeleton reorganization: A mathematical modeling and experimental study,” Journal of Theoretical Biology, vol. 202, no. 1, pp. 33–41, 2000. View at: Publisher Site  Google Scholar
 J. Pourati, A. Maniotis, D. Spiegel et al., “Is cytoskeletal tension a major determinant of cell deformability in adherent endothelial cells?” American Journal of PhysiologyCell Physiology, vol. 274, no. 5, pp. C1283–C1289, 1998. View at: Publisher Site  Google Scholar
 N. Wang, I. M. ToliNørrelykke, J. Chen et al., “Cell prestress. I. Stiffness and prestress are closely associated in adherent contractile cells,” American Journal of PhysiologyCell Physiology, vol. 282, no. 3, pp. 606–616, 2002. View at: Publisher Site  Google Scholar
 D. Stamenović and D. E. Ingber, “Models of cytoskeletal mechanics of adherent cells,” Biomechanics and Modeling in Mechanobiology, vol. 1, no. 1, pp. 95–108, 2002. View at: Publisher Site  Google Scholar
 D. E. Discher, D. H. Boal, and S. K. Boey, “Simulations of the erythrocyte cytoskeleton at large deformation. II. Micropipette aspiration,” Biophysical Journal, vol. 75, no. 3, pp. 1584–1597, 1998. View at: Publisher Site  Google Scholar
 J. Li, M. Dao, C. T. Lim, and S. Suresh, “Spectrinlevel modeling of the cytoskeleton and optical tweezers stretching of the erythrocyte,” Biophysical Journal, vol. 88, no. 5, pp. 3707–3719, 2005. View at: Publisher Site  Google Scholar
 M. A. Wozniak and C. S. Chen, “Mechanotransduction in development: a growing role for contractility,” Nature Reviews Molecular Cell Biology, vol. 10, no. 1, pp. 34–43, 2009. View at: Publisher Site  Google Scholar
 X. J. Fan, “Cell biomechanics,” Advances in Mechanics, vol. 25, no. 2, pp. 197–208, 1995. View at: Google Scholar
 R. M. Hochmuth, “Micropipette aspiration of living cells,” Journal of Biomechanics, vol. 33, no. 1, pp. 15–22, 2000. View at: Publisher Site  Google Scholar
 B. Fabry, G. N. Maksym, J. P. Butler, M. Glogauer, D. Navajas, and J. J. Fredberg, “Scaling the microrheology of living cells,” Physical Review Letters, vol. 87, no. 14, Article ID 148102, 2001. View at: Google Scholar
 V. C. Mow, S. C. Kuei, W. M. Lai, and C. G. Armstrong, “Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments,” Journal of Biomechanical Engineering, vol. 102, no. 1, pp. 73–84, 1980. View at: Publisher Site  Google Scholar
 F. P. T. Baaijens, W. R. Trickey, T. A. Laursen, and F. Guilak, “Large deformation finite element analysis of micropipette aspiration to determine the mechanical properties of the chondrocyte,” Annals of Biomedical Engineering, vol. 33, no. 4, pp. 494–501, 2005. View at: Publisher Site  Google Scholar
 J. Z. Wu and W. Herzog, “Analysis of the mechanical behavior of chondrocytes in unconfined compression tests for cyclic loading,” Journal of Biomechanics, vol. 39, no. 4, pp. 603–616, 2006. View at: Publisher Site  Google Scholar
 E. J. Koay, A. C. Shieh, and K. A. Athanasiou, “Creep indentation of single cells,” Journal of Biomechanical Engineering, vol. 125, no. 3, pp. 334–341, 2003. View at: Publisher Site  Google Scholar
 E. M. Darling, S. Zauscher, and F. Guilak, “Viscoelastic properties of zonal articular chondrocytes measured by atomic force microscopy,” Osteoarthritis and Cartilage, vol. 14, no. 6, pp. 571–579, 2006. View at: Publisher Site  Google Scholar
 V. Lulevich, T. Zink, H.Y. Chen, F.T. Liu, and G.Y. Liu, “Cell mechanics using atomic force microscopybased singlecell compression,” Langmuir, vol. 22, no. 19, pp. 8151–8155, 2006. View at: Publisher Site  Google Scholar
 A. R. Bausch, F. Ziemann, A. A. Boulbitch, K. Jacobson, and E. Sackmann, “Local measurements of viscoelastic parameters of adherent cell surfaces by magnetic bead microrheometry,” Biophysical Journal, vol. 75, no. 4, pp. 2038–2049, 1998. View at: Publisher Site  Google Scholar
 J. Guck, S. Schinkinger, B. Lincoln et al., “Optical deformability as an inherent cell marker for testing malignant transformation and metastatic competence,” Biophysical Journal, vol. 88, no. 5, pp. 3689–3698, 2005. View at: Publisher Site  Google Scholar
 J. Cao, S. Usami, and C. Dong, “Development of a sideview chamber for studying cellsurface adhesion under flow conditions,” Annals of Biomedical Engineering, vol. 25, no. 3, pp. 573–580, 1997. View at: Publisher Site  Google Scholar
 C. R. Jacobs, S. Temiyasathit, and A. B. Castillo, “Osteocyte mechanobiology and pericellular mechanics,” Annual Review of Biomedical Engineering, vol. 12, pp. 369–400, 2010. View at: Publisher Site  Google Scholar
 L. F. Bonewald and M. L. Johnson, “Osteocytes, mechanosensing and Wnt signaling,” Bone, vol. 42, no. 4, pp. 606–615, 2008. View at: Publisher Site  Google Scholar
 L. You, S. Temiyasathit, P. Lee et al., “Osteocytes as mechanosensors in the inhibition of bone resorption due to mechanical loading,” Bone, vol. 42, no. 1, pp. 172–179, 2008. View at: Publisher Site  Google Scholar
 S. C. Manolagas, “Birth and death of bone cells: basic regulatory mechanisms and implications for the pathogenesis and treatment of osteoporosis,” Endocrine Reviews, vol. 21, no. 2, pp. 115–137, 2000. View at: Publisher Site  Google Scholar
 A. D. Baik, X. L. Lu, J. Qiu et al., “Quasi3D cytoskeletal dynamics of osteocytes under fluid flow,” Biophysical Journal, vol. 99, no. 9, pp. 2812–2820, 2010. View at: Publisher Site  Google Scholar
 J. Qiu and F.F. Li, “Mechanical behavior of an individual adherent MLOY4 osteocyte under shear flow,” Biomechanics and Modeling in Mechanobiology, vol. 16, no. 1, pp. 1–12, 2017. View at: Publisher Site  Google Scholar
 S. Yamada, D. Wirtz, and S. C. Kuo, “Mechanics of living cells measured by laser tracking microrheology,” Biophysical Journal, vol. 78, no. 4, pp. 1736–1747, 2000. View at: Publisher Site  Google Scholar
 E. Evans, N. Mohandas, and A. Leung, “Static and dynamic rigidities of normal and sickle erythrocytes. Major influence of cell hemoglobin concentration,” The Journal of Clinical Investigation, vol. 73, no. 2, pp. 477–488, 1984. View at: Publisher Site  Google Scholar
 I. M. Ward and P. R. Pinnock, “The mechanical properties of solid polymers,” British Journal of Applied Physics, vol. 17, no. 1, pp. 3–32, 1966. View at: Publisher Site  Google Scholar
 W. R. Trickey, G. M. Lee, and F. Guilak, “Viscoelastic properties of chondrocytes from normal and osteoarthritic human cartilage,” Journal of Orthopaedic Research, vol. 18, no. 6, pp. 891–898, 2000. View at: Publisher Site  Google Scholar
 M. Sato, D. P. Theret, L. T. Wheeler, N. Ohshima, and R. M. Nerem, “Application of the micropipette technique to the measurement of cultured porcine aortic endothelial cell viscoelastic properties,” Journal of Biomechanical Engineering, vol. 112, no. 3, pp. 263–268, 1990. View at: Publisher Site  Google Scholar
 S. Kuroshima, M. Kaku, T. Ishimoto, M. Sasaki, T. Nakano, and T. Sawase, “A paradigm shift for bone quality in dentistry: A literature review,” Journal of Prosthodontic Research, vol. 61, no. 4, pp. 353–362, 2017. View at: Publisher Site  Google Scholar
 R. F. M. van Oers, H. Wang, and R. G. Bacabac, “Osteocyte shape and mechanical loading,” Current Osteoporosis Reports, vol. 13, no. 2, pp. 61–66, 2015. View at: Publisher Site  Google Scholar
 R. Y. Kwon and C. R. Jacobs, “Timedependent deformations in bone cells exposed to fluid flow in vitro: investigating the role of cellular deformation in fluid flowinduced signaling,” Journal of Biomechanics, vol. 40, no. 14, pp. 3162–3168, 2007. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Jun Qiu and FangFang Li. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.