Mathematical Problems in Engineering
Volume 2019 (2019), Article ID 5868561, 13 pages
https://doi.org/10.1155/2019/5868561

Research Article

Identification of Viscoelastic Constitutive Parameters of a Cell Based on Fluid-Structure Coupled Finite Element Model and Experiment

Jun Qiu1,2 and Fang-Fang Li³

¹State Key Laboratory of Hydroscience & Engineering, Tsinghua University, Beijing 100084, China
²State Key Laboratory of Plateau Ecology and Agriculture, Qinghai University, Xining 810016, China
³College of Water Resources & Civil Engineering, China Agricultural University, Beijing 100083, China

Correspondence should be addressed to Fang-Fang Li; liff@cau.edu.cn

Received 20 December 2018; Revised 11 March 2019; Accepted 26 March 2019; Published 11 April 2019

Guest Editor: Juan Melchor

Copyright © 2019 Jun Qiu and Fang-Fang 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.

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 fluid-structure 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 in-depth 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 power-law, 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], cyto-indentation [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 bone-forming osteoblast and bone-removing 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 fluid-structure 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)

Figure 1: Cell flow chamber experiment. (a) Schematic diagram of flow chamber experiment; the side images of the osteocyte (b) before and (c) after the shear flow.

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 self-built quasi-three-dimensional 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 Fluid-Structure 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 fluid-structure interaction problems. In the process of fluid-solid 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 fluid-structure 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 fluid-structure coupling solver (ADINA-FSI) 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 fluid-solid 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 cross-section tube is the basis of the interaction of cell shear flow. The numerical simulation of the flow in a square cross-section 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 cross-section height of the square tube is H and the width is W. For steady laminar flow of incompressible fluid, the Navier-Stokes equation under unidirectional flow in a square cross-section 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.

Table 1: Convergence analysis of 2D fluid meshes density.

Analytical solution
Coarse meshes
General meshes
Fine meshes

Numerical solution
Error
Numerical solution
Error
Numerical solution
Error

Max shear stress (Pa)
1.9643E-01
1.8990E-01
-3.32%
1.9320E-01
-1.64%
1.9480E-01
-0.83%

Max flow rate
(mm/s)
2.0776E-02
2.0780E-02
0.02%
2.0780E-02
0.02%
2.0780E-02
0.02%

volume metric
(mm³/s)
7.6179E-06
7.6094E-06
-0.11%
7.6158E-06
-0.03%
7.6174E-06
-0.01%

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 Fluid-Structure Coupled Finite Element Model

The fluid-structure 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.

Figure 2: 2D fluid-structure coupled finite element model under share flow.

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 Fluid-Structure Coupled Finite Element Model

2.3.1. Fluid Model

In the finite element simulation, the establishment of such full-scale 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 cross-section 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 fluid-structure coupling boundary condition is applied to the interface between fluid and cell.

2.3.2. Structure Model

Before establishing the final fluid-structure 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 fluid-structure coupling model of hemispherical cell configuration, the fluid domain is shown in Figure 3(a). The length of the square cross-section 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², 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.

Figure 3: 3D fluid-structure coupled finite element model under share flow. (a) Fluid field of the fluid-structure coupled finite element model based on hemispherical cell; (b) meshing scheme of the fluid field near the cell; (c) standard linear solid (SLS) model.

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 stress-strain relationship in time domain and applying the boundary conditions under stress relaxation conditions, the stress-strain 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.

Figure 4: Comparison of theoretical and numerical solutions for stress relaxation problem.

2.5. Identification of Viscoelastic Parameters of a Cell

A fluid-structure 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 fluid-structure 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 fluid-structure coupled finite element model. Meanwhile, the deformation process of cells under the fluid is observed by a self-built quasi-3D 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.

Figure 5: Framework of methodology to determine viscoelastic parameters of a cell.

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 Fluid-Structure Coupled Finite Element Model

2.6.1. Calibration of 2D Model

Table 2 shows the calculated values of the fluid-structure 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 first-order 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.

Table 2: Convergence analysis of fluid-structure interaction model of cell shear flow with three meshing methods.

0.4μm
0.8μm
1.6μm

Calculated Value
Calculated Value
Error
Calculated Value
Error

Horizontal Displacement (μm)
0.0772
0.0753
-2.46%
0.0711
-5.58%

Vertical Displacement (μm)
-0.0121
-0.012
-0.83%
-0.0119
-0.83%

Shear Stress (Pa)
2.2
2.03
-7.73%
1.77
-12.81%

Maximal Stress (Pa)
3.59
3.7
3.06%
3.86
4.32%

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)

Figure 6: Convergence analysis of 5 different mesh sizes. (a) Pressure distribution on the cell surface; (b) shear stress distribution on the cell surface.

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 fluid-structure 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.

Table 3: Convergence analysis of fluid-structure interaction model of cell shear flow with 5 different mesh sizes.

0.1μm
0.2μm
0.4μm
0.8μm
1.6μm

Calculated Value
Calculated Value
Error
Calculated Value
Error
Calculated Value
Error
Calculated Value
Error

Shear Stress (Pa)
2.38
2.3
-3.36%
2.2
-4.35%
2.03
-7.73%
1.77
-12.81%

Maximal Pressure (Pa)
3.56
3.58
0.56%
3.59
0.28%
3.7
3.06%
3.86
4.32%

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 fluid-structure 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.

Figure 7: The deformation of the cell under shear stress in experiments (a) and in the fluid-structure coupling finite element model with different views (b~g).

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)

Figure 8: The influence of viscoelastic parameters on the cell mechanical response under fluid shear stress: (a) is changed; (b) is changed; (c) is changed.

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)

Figure 9: Comparison between the displace response of the cell in the fluid-structure coupled finite element model and those from the experimental results. (a) The displacement at the top area of the cell along the flow direction from the experimental results for #01 osteocyte; (b) using first group of viscoelastic parameters M1 for #01 osteocyte; (c) the process of three optimization iterations; M1, M2, and M3 are three groups of decision variables for #01 osteocyte; (d) comparison results for three different osteocytes (#01, #02, and #03).

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 fluid-structure 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.

Table 4: Identification results of three viscoelastic parameters for osteocytes.

Cell#01
Cell#02
Cell#03
Average

(n=3)

(kPa)
0.11
0.17
0.18
0.15±0.02

(kPa)
0.70
0.35
0.80
0.62±0.14

(kPa⋅s)
0.85
1.3
2.0
1.38±0.33

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 fluid-solid 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 quasi-3D microscopy system. Meanwhile, a fluid-structure 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 fluid-structure 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 real-time 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.

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.

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.

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.

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.

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.

J. Pourati, A. Maniotis, D. Spiegel et al., “Is cytoskeletal tension a major determinant of cell deformability in adherent endothelial cells?” American Journal of Physiology-Cell Physiology, vol. 274, no. 5, pp. C1283–C1289, 1998.

N. Wang, I. M. Toli-Nørrelykke, J. Chen et al., “Cell prestress. I. Stiffness and prestress are closely associated in adherent contractile cells,” American Journal of Physiology-Cell Physiology, vol. 282, no. 3, pp. 606–616, 2002.

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.

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.

J. Li, M. Dao, C. T. Lim, and S. Suresh, “Spectrin-level modeling of the cytoskeleton and optical tweezers stretching of the erythrocyte,” Biophysical Journal, vol. 88, no. 5, pp. 3707–3719, 2005.

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.

X. J. Fan, “Cell biomechanics,” Advances in Mechanics, vol. 25, no. 2, pp. 197–208, 1995.

R. M. Hochmuth, “Micropipette aspiration of living cells,” Journal of Biomechanics, vol. 33, no. 1, pp. 15–22, 2000.

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.

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.

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.

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.

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.

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.

V. Lulevich, T. Zink, H.-Y. Chen, F.-T. Liu, and G.-Y. Liu, “Cell mechanics using atomic force microscopy-based single-cell compression,” Langmuir, vol. 22, no. 19, pp. 8151–8155, 2006.

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.

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.

J. Cao, S. Usami, and C. Dong, “Development of a side-view chamber for studying cell-surface adhesion under flow conditions,” Annals of Biomedical Engineering, vol. 25, no. 3, pp. 573–580, 1997.

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.

L. F. Bonewald and M. L. Johnson, “Osteocytes, mechanosensing and Wnt signaling,” Bone, vol. 42, no. 4, pp. 606–615, 2008.

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.

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.

A. D. Baik, X. L. Lu, J. Qiu et al., “Quasi-3D cytoskeletal dynamics of osteocytes under fluid flow,” Biophysical Journal, vol. 99, no. 9, pp. 2812–2820, 2010.

J. Qiu and F.-F. Li, “Mechanical behavior of an individual adherent MLO-Y4 osteocyte under shear flow,” Biomechanics and Modeling in Mechanobiology, vol. 16, no. 1, pp. 1–12, 2017.

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.

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.

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.

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.

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.

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.

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.

R. Y. Kwon and C. R. Jacobs, “Time-dependent deformations in bone cells exposed to fluid flow in vitro: investigating the role of cellular deformation in fluid flow-induced signaling,” Journal of Biomechanics, vol. 40, no. 14, pp. 3162–3168, 2007.