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

Authors: Jun Qiu; Fang-Fang Li
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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 bothk
1 and k
2 of an osteocyte are in the order of several hundred Pa and η is in the order of kPa·s.


---

## Body

## 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 bothk
1 and k
2 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 Figure1(a).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.
(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/c
m
2. 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.2D 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.2D 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](1)
0
=
-
d
P
d
x
+
μ
∂
2
u
x
∂
y
2
+
∂
2
u
x
∂
z
2where μ is the fluid viscous coefficient, P is the fluid pressure, u is the flow velocity, x is the dimension of the flow field direction, y is the direction of the height of the flow chamber, and z is the direction of the width of the flow chamber. The boundary conditions for this problem are shown in (2)
u
x
=
0
y
=
±
H
2
u
x
=
0
z
=
±
w
2where H is the height of the flow chamber and w 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(3)
u
x
=
1
2
μ
-
d
P
d
x
×
H
2
4
-
z
2
-
8
H
2
π
3
×
∑
n
=
0
∞
-
1
n
cos
⁡
2
n
+
1
π
z
/
H
cosh
⁡
2
n
+
1
π
y
/
H
2
n
+
1
3
cosh
⁡
2
n
+
1
π
w
/
2
HThe wall shear stresst
w is as in (4)
τ
w
=
-
μ
d
u
x
d
z
z
=
-
H
/
2
=
-
H
2
-
d
P
d
x
×
1
-
8
π
2
∑
n
=
0
∞
c
o
s
h
2
n
+
1
π
y
/
H
2
n
+
1
2
c
o
s
h
2
n
+
1
π
w
/
2
HThe flow rateQ in the square section is as in (5)
Q
=
w
H
3
4
μ
-
d
P
d
x
×
1
3
-
64
H
π
5
w
∑
n
=
0
∞
tanh
⁡
2
n
+
1
π
w
/
2
H
2
n
+
1
5In 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 P
a, 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.Table1 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(mm3/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 Figure2. 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.3D 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 numberR
e is as in (6)
R
e
=
ρ
v
L
μIn the cell flow chamber experiment, the flow velocityv 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 L, 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 L
ε required for the full development of the flow field, as shown in (7)
L
ε
d
=
0.59
+
0.055
R
eIn 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 Figure3(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.06mm/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 10dyn/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.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(8)
σ
t
ε
0
=
E
r
e
l
t
=
k
1
+
k
2
e
-
t
/
τ
,
τ
=
η
k
2where E
r
e
l is the relaxation modulus of the viscoelastic materials, σ
(
t
) is creep stress, and ε
0 is initial strain. It can be seen from the expression of relaxation modulus that, at the initial moment t=0, the elastic modulus is k
1
+
k
2, 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 k
1, which is called the equilibrium modulus. The viscous coefficient and k
2 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(9)
s
i
j
t
=
2
G
0
e
i
j
t
+
2
∫
0
t
e
i
j
t
-
τ
d
G
τ
d
τ
d
τ
(10)
σ
k
k
t
=
3
K
0
ε
k
k
t
+
3
∫
0
t
ε
k
k
t
-
τ
d
K
τ
d
τ
d
τ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 (11)
G
t
=
G
∞
+
∑
i
=
1
η
G
G
i
e
-
β
i
t
(12)
K
t
=
K
∞
+
∑
i
=
1
η
K
K
i
e
-
γ
i
tThe 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(13)
G
∞
=
1
3
k
1
G
1
=
1
3
k
2
β
1
=
1
τ
=
k
2
η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 (14)
σ
t
=
ε
0
k
e
+
k
1
e
-
t
/
τ
,
τ
=
η
k
1According 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(
k
1
,
k
2
,
η
), 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 Figure5. 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

Table2 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 Figure6(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.Figure 6
Convergence analysis of 5 different mesh sizes. (a) Pressure distribution on the cell surface; (b) shear stress distribution on the cell surface.
(a)
(b)As shown in Figure6(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 Table3, 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 Figures7(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. Figures7(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 k
1, the other two parameters are kept constant, k
2 = 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 k
1 significantly affects the response of the cell after equilibrium state, and k
1 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 k
1.Figure 8
The influence of viscoelastic parameters on the cell mechanical response under fluid shear stress: (a)k
1 is changed; (b) k
2 is changed; (c) η is changed.
(a)
(b)
(c)Adjusting the value ofk
2 and keeping the other two parameters constant, k
1=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 k
2 significantly affects the initial deformation response of the cell. From the second half of the response curve, it can be seen that k
2 basically does not affect the response of the cell after equilibrium state. Thus the initial deformation response can be reduced by increasing k
2.Adjusting the value of the viscosity coefficientη and keeping the other two parameters constant, k
1=120 Pa and k
2=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.


## 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 Figure1(a).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.
(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/c
m
2. 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.2D 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.2D 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](1)
0
=
-
d
P
d
x
+
μ
∂
2
u
x
∂
y
2
+
∂
2
u
x
∂
z
2where μ is the fluid viscous coefficient, P is the fluid pressure, u is the flow velocity, x is the dimension of the flow field direction, y is the direction of the height of the flow chamber, and z is the direction of the width of the flow chamber. The boundary conditions for this problem are shown in (2)
u
x
=
0
y
=
±
H
2
u
x
=
0
z
=
±
w
2where H is the height of the flow chamber and w 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(3)
u
x
=
1
2
μ
-
d
P
d
x
×
H
2
4
-
z
2
-
8
H
2
π
3
×
∑
n
=
0
∞
-
1
n
cos
⁡
2
n
+
1
π
z
/
H
cosh
⁡
2
n
+
1
π
y
/
H
2
n
+
1
3
cosh
⁡
2
n
+
1
π
w
/
2
HThe wall shear stresst
w is as in (4)
τ
w
=
-
μ
d
u
x
d
z
z
=
-
H
/
2
=
-
H
2
-
d
P
d
x
×
1
-
8
π
2
∑
n
=
0
∞
c
o
s
h
2
n
+
1
π
y
/
H
2
n
+
1
2
c
o
s
h
2
n
+
1
π
w
/
2
HThe flow rateQ in the square section is as in (5)
Q
=
w
H
3
4
μ
-
d
P
d
x
×
1
3
-
64
H
π
5
w
∑
n
=
0
∞
tanh
⁡
2
n
+
1
π
w
/
2
H
2
n
+
1
5In 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 P
a, 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.Table1 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(mm3/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 Figure2. 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.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.2D 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](1)
0
=
-
d
P
d
x
+
μ
∂
2
u
x
∂
y
2
+
∂
2
u
x
∂
z
2where μ is the fluid viscous coefficient, P is the fluid pressure, u is the flow velocity, x is the dimension of the flow field direction, y is the direction of the height of the flow chamber, and z is the direction of the width of the flow chamber. The boundary conditions for this problem are shown in (2)
u
x
=
0
y
=
±
H
2
u
x
=
0
z
=
±
w
2where H is the height of the flow chamber and w 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(3)
u
x
=
1
2
μ
-
d
P
d
x
×
H
2
4
-
z
2
-
8
H
2
π
3
×
∑
n
=
0
∞
-
1
n
cos
⁡
2
n
+
1
π
z
/
H
cosh
⁡
2
n
+
1
π
y
/
H
2
n
+
1
3
cosh
⁡
2
n
+
1
π
w
/
2
HThe wall shear stresst
w is as in (4)
τ
w
=
-
μ
d
u
x
d
z
z
=
-
H
/
2
=
-
H
2
-
d
P
d
x
×
1
-
8
π
2
∑
n
=
0
∞
c
o
s
h
2
n
+
1
π
y
/
H
2
n
+
1
2
c
o
s
h
2
n
+
1
π
w
/
2
HThe flow rateQ in the square section is as in (5)
Q
=
w
H
3
4
μ
-
d
P
d
x
×
1
3
-
64
H
π
5
w
∑
n
=
0
∞
tanh
⁡
2
n
+
1
π
w
/
2
H
2
n
+
1
5In 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 P
a, 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.Table1 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(mm3/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 Figure2. 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.3D 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 numberR
e is as in (6)
R
e
=
ρ
v
L
μIn the cell flow chamber experiment, the flow velocityv 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 L, 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 L
ε required for the full development of the flow field, as shown in (7)
L
ε
d
=
0.59
+
0.055
R
eIn 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 Figure3(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.06mm/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 10dyn/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.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.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 numberR
e is as in (6)
R
e
=
ρ
v
L
μIn the cell flow chamber experiment, the flow velocityv 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 L, 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 L
ε required for the full development of the flow field, as shown in (7)
L
ε
d
=
0.59
+
0.055
R
eIn 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 Figure3(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.06mm/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 10dyn/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.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(8)
σ
t
ε
0
=
E
r
e
l
t
=
k
1
+
k
2
e
-
t
/
τ
,
τ
=
η
k
2where E
r
e
l is the relaxation modulus of the viscoelastic materials, σ
(
t
) is creep stress, and ε
0 is initial strain. It can be seen from the expression of relaxation modulus that, at the initial moment t=0, the elastic modulus is k
1
+
k
2, 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 k
1, which is called the equilibrium modulus. The viscous coefficient and k
2 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(9)
s
i
j
t
=
2
G
0
e
i
j
t
+
2
∫
0
t
e
i
j
t
-
τ
d
G
τ
d
τ
d
τ
(10)
σ
k
k
t
=
3
K
0
ε
k
k
t
+
3
∫
0
t
ε
k
k
t
-
τ
d
K
τ
d
τ
d
τ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 (11)
G
t
=
G
∞
+
∑
i
=
1
η
G
G
i
e
-
β
i
t
(12)
K
t
=
K
∞
+
∑
i
=
1
η
K
K
i
e
-
γ
i
tThe 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(13)
G
∞
=
1
3
k
1
G
1
=
1
3
k
2
β
1
=
1
τ
=
k
2
η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 (14)
σ
t
=
ε
0
k
e
+
k
1
e
-
t
/
τ
,
τ
=
η
k
1According 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(
k
1
,
k
2
,
η
), 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 Figure5. 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

Table2 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 Figure6(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.Figure 6
Convergence analysis of 5 different mesh sizes. (a) Pressure distribution on the cell surface; (b) shear stress distribution on the cell surface.
(a)
(b)As shown in Figure6(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 Table3, 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 Figures7(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. Figures7(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 k
1, the other two parameters are kept constant, k
2 = 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 k
1 significantly affects the response of the cell after equilibrium state, and k
1 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 k
1.Figure 8
The influence of viscoelastic parameters on the cell mechanical response under fluid shear stress: (a)k
1 is changed; (b) k
2 is changed; (c) η is changed.
(a)
(b)
(c)Adjusting the value ofk
2 and keeping the other two parameters constant, k
1=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 k
2 significantly affects the initial deformation response of the cell. From the second half of the response curve, it can be seen that k
2 basically does not affect the response of the cell after equilibrium state. Thus the initial deformation response can be reduced by increasing k
2.Adjusting the value of the viscosity coefficientη and keeping the other two parameters constant, k
1=120 Pa and k
2=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.


## 2.6.1. Calibration of 2D Model

Table2 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 Figure6(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.Figure 6
Convergence analysis of 5 different mesh sizes. (a) Pressure distribution on the cell surface; (b) shear stress distribution on the cell surface.
(a)
(b)As shown in Figure6(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 Table3, 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 Figures7(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. Figures7(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 k
1, the other two parameters are kept constant, k
2 = 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 k
1 significantly affects the response of the cell after equilibrium state, and k
1 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 k
1.Figure 8
The influence of viscoelastic parameters on the cell mechanical response under fluid shear stress: (a)k
1 is changed; (b) k
2 is changed; (c) η is changed.
(a)
(b)
(c)Adjusting the value ofk
2 and keeping the other two parameters constant, k
1=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 k
2 significantly affects the initial deformation response of the cell. From the second half of the response curve, it can be seen that k
2 basically does not affect the response of the cell after equilibrium state. Thus the initial deformation response can be reduced by increasing k
2.Adjusting the value of the viscosity coefficientη and keeping the other two parameters constant, k
1=120 Pa and k
2=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 Figure9(a). The displacement response is typical as viscoelastic materials, which is also similar to the deformation response curve observed by other researchers [38].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).
(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 fluid-structure coupled finite element model is shown in Figure9(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 ofk
1 and k
2 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. Figure9(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 Figure9(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 k
1 and k
2 are in the order of several hundred Pa andηis in the order ofkPa·s.Table 4
Identification results of three viscoelastic parameters for osteocytes.
Cell#01
Cell#02
Cell#03
Average
(n=3)
k
1 (kPa)
0.11
0.17
0.18
0.15±0.02
k
2 (kPa)
0.70
0.35
0.80
0.62±0.14
η (kPa⋅s)
0.85
1.3
2.0
1.38±0.33Our 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 bothk
1 and k
2 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.



---

*Source: 5868561-2019-04-11.xml*