A hyperelastic model for simulating cells in flow
Sebastian J. Müller, Franziska Weigl, Carina Bezold, Ana Sancho, Christian Bächer, Krystyna Albrecht, Stephan Gekle
PPublished 2020 in: Biomechanics and Modeling in MechanobiologyDOI: 10.1007/s10237-020-01397-2
Unfortunately, the journal version misses an important contributor.The correct author list is the one given in this document.
A hyperelastic model for simulating cells in flow
Sebastian J. M¨uller , Franziska Weigl , Carina Bezold , Ana Sancho , ,Christian B¨acher , Krystyna Albrecht and Stephan Gekle Theoretical Physics VI, Biofluid Simulation and Modeling, University of Bayreuth,Universit¨atsstraße 30, 95440 Bayreuth, Germany Department of Functional Materials in Medicine and Dentistry and Bavarian PolymerInstitute (BPI), University of W¨urzburg, Pleicherwall 2, 97070 W¨urzburg, Germany Department of Automatic Control and Systems Engineering, University of the BasqueCountry UPV/EHU, San Sebastian, Spain
Abstract
In the emerging field of 3D bioprinting, cell damage due to large deformations is considered a main cause for celldeath and loss of functionality inside the printed construct. Those deformations, in turn, strongly depend on themechano-elastic response of the cell to the hydrodynamic stresses experienced during printing. In this work, wepresent a numerical model to simulate the deformation of biological cells in arbitrary three-dimensional flows.We consider cells as an elastic continuum according to the hyperelastic Mooney–Rivlin model. We then employforce calculations on a tetrahedralized volume mesh.To calibrate our model, we perform a series of FluidFM ® compression experiments with REF52 cells demon-strating that all three parameters of the Mooney–Rivlin model are required for a good description of the experi-mental data at very large deformations up to 80 %. In addition, we validate the model by comparing to previousAFM experiments on bovine endothelial cells and artificial hydrogel particles. To investigate cell deformationin flow, we incorporate our model into Lattice Boltzmann simulations via an Immersed-Boundary algorithm. Inlinear shear flows, our model shows excellent agreement with analytical calculations and previous simulationdata. Keywords: Hyperelasticity, Cell deformation, Mooney–Rivlin, Atomic force Microscopy, Shear flow, Lattice-Boltzmann
The dynamic behavior of flowing cells is central to thefunctioning of organisms and forms the base for a va-riety of biomedical applications. Technological sys-tems that make use of the elastic behavior of cells are,for example, cell sorting [1], real-time deformabilitycytometry [2, 3] or probing techniques for cytoskele-tal mechanics [4–15]. In most, but not all, of theseapplications cell deformations typically remain rathersmall. A specific example where large deformationsbecome important is 3D bioprinting. Bioprinting is atechnology which, analogously to common 3D print-ing, pushes a suspension of cells in highly viscous hydrogels—a so-called bioink—through a fine nozzleto create three-dimensional tissue structures. A majorchallenge in this process lies in the control of large celldeformations and cell damage during printing. Thosedeformations arise from hydrodynamic stresses in theprinter nozzle and ultimately affect the viability andfunctionality of the cells in the printed construct [16–20]. How exactly these hydrodynamic forces corre-late with cell deformation, however, strongly dependson the elastic behavior of the cell and its interactionwith the flowing liquid. Theoretical and computa-tional modeling efforts in this area have thus far beenrestricted to pure fluid simulations without actually in- a r X i v : . [ phy s i c s . b i o - ph ] D ec THEORY corporating the cells [17, 21, 22] or simple 2D geome-tries [23, 24]. The complexity of cell mechanics andthe diversity of possible applications make theoreticalmodeling of cell mechanics in flow a challenge which,to start with, requires reliable experimental data forlarge cell deformations.The most appropriate tool to measure cellular re-sponse at large deformations is atomic force mi-croscopy (AFM) [8, 25–34]. AFM cantilevers withpyramidal tips, colloidal probes, or flat geometriesare used to indent or compress cells. Therefore, acommon approach to characterize the elasticity ofcells utilizes the Hertzian theory, which describes thecontact between two linear elastic solids [35, p. 90-104], but is limited to the range of small deforma-tions [36]. Experimental measurements with medium-to-large deformations typically show significant de-viations from the Hertz prediction, e. g., for cells orhydrogel particles [37]. Instead of linear elasticity,a suitable description of cell mechanics for bioprint-ing applications requires more advanced hyperelas-tic material properties. While for simple anucleatefluid-filled cells such as, e.g., red blood cells, the-oretical models abound [38–42], the availability ofmodels for cells including a complex cytoskeleton israther limited. In axisymmetric geometries, Cailleet al. [43] and Mokbel et al. [44] used an axisymmetricfinite element model with neo-Hookean hyperelastic-ity to model AFM and microchannel experiments onbiological cells. In shear flow, approximate analyti-cal treatments are possible [45–48]. Computationally,Gao and Hu [46] carried out 2D simulations while in3D Lykov et al. [49] utilized a DPD technique basedon a bead-spring model. Furthermore, Villone et al.[50, 51] presented an arbitrary Lagrangian-Eulerianapproach for elastic particles in viscoelastic fluids. Fi-nally, Rosti et al. [52] and Saadat et al. [53] consideredviscoelastic and neo-Hookean finite element models,respectively, in shear flow.In this work, we introduce and calibrate a com-putational model for fully three-dimensional simula-tions of cells in arbitrary flows. Our approach usesa Lattice-Boltzmann solver for the fluid and a directforce formulation for the elastic equations. In contrastto earlier works [43, 44, 47, 52, 53] our model usesa three-parameter Mooney–Rivlin elastic energy func-tional. To demonstrate the need for this more com-plex elastic model, we carry out extensive FluidFM ® indentation experiments for REF52 (rat embryonic fi-broblast) cells at large cell deformation up to 80 %[54]. In addition, our model compares favorably with previous AFM experiments on bovine endothe-lial cells [43] as well as artificial hydrogel particles[37]. Our model provides a much more realisticforce–deformation behavior compared to the small-deformation Hertz approximation, but is still simpleand fast enough to allow the simulation of dense cellsuspensions in reasonable time. Particularly, our ap-proach is less computationally demanding than con-ventional finite-element methods which usually re-quire large matrix operations. Furthermore, it is eas-ily extensible and allows, e.g., the inclusion of a cellnucleus by the choice of different elastic moduli fordifferent parts of the volume.We finally present simulations of our cell model indifferent flow scenarios using an Immersed-Boundaryalgorithm to couple our model with Lattice Boltzmannfluid calculations. In a plane Couette (linear shear)flow, we investigate the shear stress dependency ofsingle cell deformation, which we compare to the av-erage cell deformation in suspensions with higher vol-ume fractions, and show that our results in the neo-Hookean limit are in accordance with earlier elasticcell models [47, 52, 53]. In general, hyperelastic models are used to describematerials that respond elastically to large deforma-tions [55, p. 93]. Many cell types can be subjected tolarge reversible shape changes. This section providesa brief overview of the hyperelastic Mooney–Rivlinmodel implemented in this work.The displacement of a point is given by u i = y i − x i , (1)where x i ( i = , ,
3) refers to the undeformed config-uration (material frame) and y i to the deformed co-ordinates (spatial frame). We define the deformationgradient tensor and its inverse as [55, p. 14,18] F i j = ∂ y i ∂ x j = ∂ u i ∂ x j + δ i j and F − i j = ∂ x i ∂ y j . (2)Together with the right Cauchy-Green deformationtensor, C = F (cid:124) F (material description), we can de-fine the following invariants which are needed for thestrain energy density calculation below: J = det F (3) I = T C J − / (4) K = (cid:0) T C − T C (cid:1) J − / (5)2 TETRAHEDRALIZED CELL MODEL
Here, T C = tr C and T C = tr (cid:0) C (cid:1) (6)are the trace of the right Cauchy-Green deformationtensor and its square, respectively. The nonlinearstrain energy density of the Mooney–Rivlin model isgiven by [56, 57] U = (cid:104) µ ( I − ) + µ ( K − ) + κ ( J − ) (cid:105) , (7)where µ , µ , and κ are material properties. Theycorrespond—for consistency with linear elasticity inthe range of small deformations—to the shear modu-lus µ = µ + µ and bulk modulus κ of the materialand are therefore related to the Young’s modulus E and the Poisson ratio ν via [55, p. 74] µ = E ( + ν ) and κ = E ( − ν ) . (8)Through the choice µ = U NH = (cid:104) µ ( I − ) + κ ( J − ) (cid:105) (9)As we show later, this can be a sufficient descriptionfor some cell types. To control the strength of the sec-ond term and quickly switch between neo-Hookeanand Mooney–Rivlin strain energy density calculation,we introduce a factor w ∈ [ , ] and set µ = w µ and µ = ( − w ) µ (10)such that w =
1, which equals setting µ = w < µ -termand thus leads to a more pronounced strain hardeningas shown in figure S-6 of the Supporting Information. In this section we apply the hyperelastic theory of sec-tion 2 to a tetrahedralized mesh as shown in figure 1.
We consider a mesh consisting of tetrahedral elementsas depicted in figure 1. The superscript α refers tothe four vertices of the tetrahedron. The elastic forceacting on vertex α in direction i is obtained from (7) by differentiating the strain energy density U with respectto the vertex displacement as f α i = − V ∂ U ∂ u α i , (11)where V is the reference volume of the tetrahedron.In contrast to Saadat et al. [53], the numerical calcu-lation of the force in our model does not rely on theintegration of the stress tensor, but on a differentiationwhere the calculation of all resulting terms involvesonly simple arithmetics. Applying the chain rule fordifferentiation yields: f α i = − V (cid:34) (cid:18) ∂ U ∂ I ∂ I ∂ T C + ∂ U ∂ K ∂ K ∂ T C (cid:19) ∂ T C ∂ F kl + (cid:18) ∂ U ∂ I ∂ I ∂ J + ∂ U ∂ K ∂ K ∂ J + ∂ U ∂ J (cid:19) ∂ J ∂ F kl + ∂ U ∂ K ∂ K ∂ T C ∂ T C ∂ F kl (cid:35) ∂ F kl ∂ u α i (12)The evaluation of (12) requires the calculation of thedeformation gradient tensor F , which is achieved bylinear interpolation of the coordinates and displace-ments inside each tetrahedral mesh element as detailedin the next section. We note that our elastic force cal-culation is purely local making it straightforward toemploy different elastic models in different regions ofthe cell and/or to combine it with elastic shell models.This flexibility can be used to describe, e.g., the cellnucleus [43] or an actin cortex [58] surrounding thecell interior. Following standard methods, e.g. Bower [55], we startby interpolating a point x i inside a single tetrahedronusing the vertex positions x α i ( α = , , , ( ξ , ξ , ξ ) with 0 ≤ ξ i ≤ , asdepicted in figure 1a. One vertex defines the originwhile the remaining three indicate the coordinate axes.A set of shape functions, i. e., interpolation functions, N α ( ξ , ξ , ξ ) is employed to interpolate positions in-side the tetrahedron volume. An arbitrary point x i in-side the element is interpolated as x i = ∑ α = N α ( ξ , ξ , ξ ) x α i , (13) Bower [55, p. 481,483] erroneously states a range of − ≤ ξ i ≤ .3 Taylor deformation parameter 3 TETRAHEDRALIZED CELL MODEL where the shape functions are defined as [55, p. 483]: N ( ξ , ξ , ξ ) = ξ (14) N ( ξ , ξ , ξ ) = ξ (15) N ( ξ , ξ , ξ ) = ξ (16) N ( ξ , ξ , ξ ) = − ξ − ξ − ξ (17)According to (1), the displacement of vertex α in i -direction is given by u α i = y α i − x α i . (18)Therefore similar to (13), the displacement at an ar-bitrary point in the volume can also be expressed interms of the shape functions and the vertex displace-ments as u i = ∑ α = N α ( ξ , ξ , ξ ) u α i . (19)The calculation of the deformation gradient tensor ac-cording to (2) requires the spatial derivative of the dis-placement: F i j − δ i j = ∂ u i ∂ x j = ∂ u i ∂ ξ k ∂ ξ k ∂ x j = A ik B k j (20)By inserting (19) into (20) and evaluating the shapefunctions, the components of the matrix A are easilydetermined to be the difference of the displacementsbetween the origin (vertex 4) and the remaining ver-tices 1, 2 and 3: A ik = u ki − u i (21)Note that due to the linear interpolation A ik is constantinside a given tetrahedron. The matrix B = J − is theinverse of the Jacobian matrix, obtained similarly to(21) as J ik = ∂ x i ∂ ξ k = x ki − x i . (22)Since x i refers to the reference coordinates, the cal-culation of the matrices J and B has to be performedonly once at the beginning of a simulation. With theinterpolation of the displacement in each tetrahedron,we can write all derivatives occurring in (12), as listedin the following: ∂ U ∂ I = µ ∂ I ∂ T C = J − ∂ U ∂ K = µ ∂ K ∂ T C = T C J − ∂ T C ∂ F il = F il ∂ I ∂ J = − T C J − ∂ K ∂ J = − (cid:0) T C − T C (cid:1) J − ∂ U ∂ J = κ ( J − ) ∂ J ∂ F il = J F − li ∂ K ∂ T C = − J − ∂ T C ∂ F il = F ik C kl ∂ F kl ∂ u α i = δ ki B ml ( δ m α − δ α ) As a measure for the cell deformation, we use the Tay-lor deformation parameter [53, 59–61] D = a − a a + a , (23)where a and a are respectively the minor and majorsemi axis of an ellipsoid corresponding to the inertiatensor of the cell. The Taylor deformation is a goodmeasure for approximately elliptic cell deformations,as they occur in shear flow (cf. section 6).To calculate D , first the components of the inertiatensor Θ i j = (cid:90) V x k x k δ i j − x i x j d V , (24)where (cid:126) x is a vector inside the volume V , are calculatedusing our discretized cell with N tet tetrahedra as Θ i j = N tet ∑ l = V l (cid:16) r lk r lk δ i j − r li r lj (cid:17) . (25)The vector (cid:126) r l denotes the center of mass of the l th tetrahedron and V l is its current volume. The eigen-values θ > θ > θ of Θ can be used to fit the semiaxes a < a < a of the corresponding ellipsoid: a = M ( − θ + θ + θ ) a = M ( θ − θ + θ ) a = M ( θ + θ − θ ) (26)The prefactor contains the mass M of the ellipsoid(considering uniform mass density) and drops out inthe calculation of D .4 COMPARISON TO FLUIDFM ® MEASUREMENTS (a) ξ ξ ξ (b) (c)Figure 1: (a) The four noded tetrahedron as mesh element within a local dimensionless coordinate system { ξ , ξ , ξ } . (b) The spherical cell model with its triangulated surface and (c) its inner tetrahedralized mesh ® measure-ments on REF52 cells In this section, we validate compression simulationsof our cell model with FluidFM ® compression exper-iments of REF52 cells stably expressing paxillin-YFP[54]. These experiments provide as an output the re-quired force to produce a certain deformation of thecell, which can be directly compared to our model.We start with a detailed description of the experimentsand show the suitability of our model to describe theelastic behavior of REF52 cells afterwards. ® indentation measure-ments We perform a series of compression measurementsof REF52 cells with a Flex FPM (Nanosurf GmbH,Germany) system that combines the AFM withthe FluidFM ® technology (Cytosurge AG, Switzer-land). In contrast to conventional AFM techniques,FluidFM ® uses flat cantilevers that possess a mi-crochannel connected to a pressure system. By ap-plying a suction pressure, cells can be aspirated andretained at the aperture of the cantilever’s tip. A moredetailed description of the setup and its functionalityis already reported in [31]. All experiments are basedon a cantilever with an aperture of 8 µm diameter anda nominal spring constant of 2 N m − . In order to mea-sure the cellular deformation, a cell was sucked ontothe tip and compressed between the cantilever and thesubstrate until a setpoint of 100 nN was reached. Im-mediately before the experiment, the cells were de- Figure 2: Example micrograph showing the FluidFM ® cantilever and a cell viewed from the top. Scale bar is30 µmtached by using Accutase (Sigma Aldrich) and weretherefore in suspension at the time of indentation. Inthis way, it can be ensured that only a single cell isdeformed during each measurement.An example micrograph of the experiment beforecompression is shown in figure 2. Analogously toAFM, primary data in form of cantilever position (inm) and deflection (in V) has to be converted to forceand deformation through the deflection sensitivity (inm V − ) and the cantilevers’ spring constant. The cel-lular deformation further requires the determination ofthe contact point, which we choose as the cantileverposition where the measured force starts to increase.The undeformed cell size is obtained as mean froma horizontal and vertical diameter measurement usingthe software imageJ.5 .2 Simulation setup 5 COMPARISON TO OTHER SETUPS The experimental setup of the previous section is eas-ily transferred and implemented for our cell model:the undeformed spherical cell rests on a fixed platewhile a second plate approaches from above to com-press the cell as depicted in figure 3 (a and b). Insection 5.2 below we will also use a slightly modi-fied version where a sphere indents the cell as shownin figure 3 (c and d). A repulsive force prevents thecell vertices from penetrating the plates or the spheri-cal indenter. The elastic restoring forces (cf. section 3)acting against this imposed compression are transmit-ted throughout the whole mesh, deforming the cell.We use meshes consisting of 2000 to 5000 verticesand about 10 000 to 30 000 tetrahedra to build up aspherical structure. More details of the mesh and itsgeneration (section S-2.4) as well as the algorithm(section S-3) are provided in the SI.
In our FluidFM ® experiment series with REF52 cells,the cell radii lie between 7 . . . ( ) µm. In figure 4 we depictthe force as function of the non-dimensionalized de-formation, i. e., the absolute compression divided bythe cell diameter. The experimental data curves sharegeneral characteristics: The force increases slowly inthe range of small deformations up to roughly 40 %,while a rapidly increasing force is observed for largerdeformations. Although the variation of the cell radiusin the different measurements is already taken into ac-count in the deformation, the point of the force upturndiffers significantly which indicates a certain variabil-ity in the elastic parameters of the individual cells.We use the compression simulation setup as de-tailed in section 4.2 to calculate force–deformationcurves of our cell model. The Poisson ratio is cho-sen as ν = .
48. In section S-2.7 of the Support-ing Information we show that variations of the ν donot strongly affect the results. A best fit approach isused to determine the Young’s modulus and the ra-tio of shear moduli w and leads to very good agree-ment between model prediction and experimental dataas shown in figure 4 as well as section S-1 of the SI.While the general range of force values is controlledusing the Young’s modulus, the Mooney–Rivlin ratio w especially defines the point of the force upturn. Wefind Young’s moduli in the range 110 Pa to 160 Pa and w = .
25, 0 .
5, and 1. For very small deformations ourhyperelastic model produces the same results as would be expected from a linear elastic model according tothe Hertz theory. See the SI (section S-2.5) for furtherdetails on the calculation of the force–deformation ac-cording to the Hertzian theory. For large deformations,the force rapidly increases due to its nonlinear charac-ter, showing strain-hardening behavior and huge devi-ations from the Hertz theory. Overall, we find an ex-cellent match between simulation and our FluidFM ® measurements with REF52 cells. In this section, we compare our simulations to ax-isymmetric calculations using the commercial soft-ware Abaqus and validate our cell model with furtherexperimental data for bovine endothelial cells from[43] and very recent data for hydrogel particles from[37].
To validate our model numerically, we compare oursimulated force–deformation curves to calculationsusing the commercial software Abaqus [62] (version6.14).In Abaqus, we use a rotationally symmetric setupconsisting of a two-dimensional semicircle, which iscompressed between two planes, similar to our simu-lation setup in section 4.2 and the finite element modelutilized in [43]. The semicircle has a radius r =
15 µm,a Young’s modulus of E = .
25 kPa and a Poisson ra-tio of ν = .
48. We choose a triangular mesh andthe built-in implementation of the hyperelastic neo-Hookean model. In figure 5 we see very good agree-ment between the results of the two different numeri-cal methods.
To compare with the AFM experiments of Caille et al.[43], we simulate a cell with radius 15 µm using thesetup of section 4.2. For the hydrogel particle inden-tation [37] we use the setup depicted in figure 3 (c andd) with a particle radius of 40 µm and a radius of thecolloidal probe of 26 . .
48 in all simulations and the Young’s modulus6
APPLICATION IN SHEAR FLOW (a) (b) (c) (d)Figure 3: (a and b) Cell compression simulations: The cell is compressed between a lower, resting, and anupper, moving, plate. (c and d) Colloidal probe cell indentation simulations: The cell rests on a plate, whilebeing indented with a sphere N o r m a l f o r ce [ n N ] Deformation [%]
Experiment Our model , w = 1120 Pa , w = 0 . , w = 0 . , w = 0 . Figure 4: Our numerical model in comparison to ourFluidFM ® measurements on REF52 cells. The labelsgive the two fit parameters E and w . We find Young’smoduli in the range of 110 Pa to 160 Pa. The Hertztheory is shown for a Young’s modulus of 180 Pa N o r m a l f o r ce [ n N ] Deformation [%]
Hertz theory3D modelAxisymmetric model
Figure 5: Comparison of force–deformation curvesobtained from our model (red line) with the linear elas-tic Hertz theory (black line) and the two-dimensionalsimulation with Abaqus (red squares), showing goodagreement between our three-dimensional and the ax-isymmetric model is determined using a best fit to the experimental datapoints. Since the neo-Hookean description appears tobe sufficient for these data sets, we further set w = ±
100 Paand the experimental data. For both systems, figure 6shows large deviations between the Hertzian theoryand the experimental data for medium-to-large defor-mations. Our model provides a significant improve-ment in this range.
We now apply our model to study the behavior of cellsin a plane Couette (linear shear) flow setup and com-pare the steady cell deformation to other numericaland analytical cell models of Gao et al. [47], Rostiet al. [52] and Saadat et al. [53]. A sketch of thesimulation setup is shown in figure 7. For simplic-ity, we choose w = µ and κ (or E and ν ), obtaining a compressible neo-Hookean form.7 .1 Single cell simulation 6 APPLICATION IN SHEAR FLOW (a) N o r m a l f o r ce [ n N ] Deformation [%]
Caille Our model12345Hertz theory2400 Pa2000 Pa1000 Pa650 Pa550 Pa (b) . .
52 0 10 20 30 40 50 60 70 80 N o r m a l f o r ce [ µ N ] Indentation [%]
Neubauer et al.
Hertz theory 469 PaOur model 580 PaOur model 580 ±
100 Pa
Figure 6: (a) Our numerical model in comparisonto experimental measurements of bovine endothelialcells from [43]. The black line depicts the predictionof the Hertz theory for a Young’s modulus of 1000 Pa.(b) Our numerical model in comparison to experi-mental measurements of hydrogel particles from [37].The indicated range corresponds to the experimentallyfound range of ±
100 Pa for the Young’s modulus ac-cording to the depicted Hertz model We use the Lattice Boltzmann implementation of theopen source software package ESPResSo [63, 64].Coupling between fluid and cell is achieved via theimmersed-boundary algorithm [53, 65] which we im-plemented into ESPResSo [58, 66]. We note here that,in contrast to Saadat et al. [53], we do not subtractthe fluid stress within the particle interior. This leadsto a small viscous response of the cell material in ad-dition to its elasticity. To obtain (approximately) thelimit of a purely elastic particle, we exploit a recentlydeveloped method by Lehmann et al. [67] to discrim-inate between the cell interior and exterior during thesimulation. Using this technique, we can tune the ra-tio between inner and outer viscosity λ with λ → λ = D (23) of our initially spherical cell model in shearflow at different shear rates ˙ γ . The first simulation setup, a single cell in infinite shearflow, is realized by choosing a simulation box of thedimensions 10 × × x × y × z ) in units of the cellradius. The infinite shear flow is approximated by ap-plying a tangential velocity u wall on the x - z -planes at y = y =
15 in positive x -direction,as depicted in figure 7. The tangential wall velocity iscalculated using the distance H of the parallel planesand the constant shear rate ˙ γ via u wall = H ˙ γ . (27)The box is periodic in x and z . A single cell is placedat the center of the simulation box corresponding to avolume fraction of φ = . ρ = kg m − ,dynamic viscosity η = − Pa s, and shear rate ˙ γ = − . The capillary number is defined by [46]Ca = η ˙ γµ , (28)and is used to set the shear modulus µ of our cell rel-ative to the fluid shear stress η ˙ γ . Simulation snap-shots of the steady state deformation of a single cell inshear flow are depicted in dependency of the capillarynumber in figure 8a. We compare the Taylor deforma-tion parameter D to previous approximate analyticalcalculations of Gao et al. [47] for a three-dimensional8 .2 Multiple cell simulations 6 APPLICATION IN SHEAR FLOW xy − u wall u wall H R
Figure 7: Schematic of the single cell in shear flow.The cell sits in the center of the box and shows anapproximately elliptic deformation as well as tank-treading, i. e., a rotation of the membrane around thesteady shape in the x - y -planeelastic solid in infinite shear flow in figure 8b and seereasonable agreement for our standard case of λ = λ = .
05, i.e.close to the limit of a purely elastic solid, the agree-ment is nearly perfect. Finally, we demonstrate thatthe elastic particle exhibits a tank-treading motion insection S-4.2.A possibly even more intuitive way to measure celldeformation is the net strain of the cell which we de-fine as ∆ ε = ( d max − d ref ) d ref . (29)It describes the relative stretching of the cell usingthe maximum elongation d max , i. e., the maximum dis-tance of two cell vertices, and its reference diameter d ref = R . A strain of ∆ ε = ∆ ε as functionof Ca. For small capillary numbers, i. e., small shearstresses, a linear stress-strain dependency is observed.Above Ca ≈ .
3, the strain-hardening, nonlinear be-havior of the neo-Hookean model can be seen. Bystretching the cell up to 280 % of its initial size, thisplot demonstrates again the capability of our model tosmoothly treat large deformations.
The second simulation setup, implemented to investi-gate the multiple particle aspect of our model, consistsof 4 (8) cells in a 5 × × φ = .
11 ( φ = .
22) occupied by cells. The cells (a) Ca = . . . . . . . . . . . . . . . . (b) . . . . . . D Ca Gao et al. λ = 1Our model λ = 0 . (c) . .
52 0 0 . . ∆ (cid:15) CaOur model λ = 1Linear fit ∆ (cid:15) (Ca) = 1 .
367 Ca
Figure 8: (a) Converged shapes of a single cell in a10 × × x × y × z ) simulation box (in units of thecell radius) with a shear flow in x -direction as func-tion of the capillary number Ca. (b) Comparison ofour model predictions for a single cell in shear flow tothe analytical 3D calculations in figure 7 of Gao et al.[47] in the range of Ca ∈ [ . , . ] . (c) The relativestretch ∆ ε of our cell model as function of the capillarynumber Ca. A linear behavior is found for small cap-illary numbers up to Ca = .
3, while increasing stressis required for larger deformations due to the strain-hardening quality of the neo-Hookean model. Linesare a guide to the eye9
CONCLUSION (a) xy z φ = . φ = . = . (b) . . . . . . . . . D Ca φ = 0 . φ = 0 . et al. Rosti et al.
Our model
Figure 9: (a) Multiple cells in a 5 × × x × y × z )simulation box (in units of the cell radius) with a con-fined shear flow in x -direction for a capillary num-ber of Ca = . φ = .
11, and 8 cells corresponding to φ = .
22. (b) Averaged deformation of multiple cellsimulations with φ = .
11 and φ = .
22 in compari-son to data from figure 3 of Rosti et al. [52] and figure13 of Saadat et al. [53]are inserted at random initial positions in the box andthe flow parameters are the same as in the first setup(cf. section 6.1).Figure 9a shows simulation snapshots of the cellsin suspensions with volume fraction φ = .
11 and φ = .
22 for Ca = .
2. The Taylor deformation ofthe suspensions, depicted in figure 9b, is calculated asan average over all cells and over time after an ini-tial transient timespan. We find good agreement whencomparing the averaged cell deformation in suspen-sion with Rosti et al. [52], Saadat et al. [53].
We presented a simple but accurate numerical modelfor cells and other microscopic particles for the use in computational fluid-particle dynamics simulations.The elastic behavior of the cells is modeled by ap-plying Mooney–Rivlin strain energy calculations ona uniformly tetrahedralized spherical mesh. We per-formed a series of FluidFM ® compression experi-ments with REF52 cells as an example for cells usedin bioprinting processes and found excellent agree-ment between our numerical model and the measure-ments if all three parameters of the Mooney–Rivlinmodel are used. In addition, we showed that the modelcompares very favorably to force versus deformationdata from previous AFM compression experiments onbovine endothelial cells [43] as well as colloidal probeAFM indentation of artificial hydrogel particles [37].At large deformations, a clear improvement comparedto Hertzian contact theory has been observed.By coupling our model to Lattice Boltzmann fluidcalculations via the Immersed-Boundary method, thecell deformation in linear shear flow as function of thecapillary number was found in good agreement withanalytical calculations by Gao et al. [47] on isolatedcells as well as previous simulations of neo-Hookeanand viscoelastic solids [52, 53] at various volume frac-tions.The presented method together with the precise de-termination of model parameters by FluidFM ® /AFMexperiments may provide an improved set of tools topredict cell deformation - and ultimately cell viability- in strong hydrodynamic flows as occurring, e.g., inbioprinting applications. Acknowledgements
Funded by the Deutsche Forschungsgemeinschaft(DFG, German Research Foundation) — Project num-ber 326998133 — TRR 225 “Biofabrication” (subpro-ject B07). We gratefully acknowledge computing timeprovided by the SuperMUC system of the LeibnizRechenzentrum, Garching. We further acknowledgesupport through the computational resources providedby the Bavarian Polymer Institute. Christian B¨acherthanks the Studienstiftung des deutschen Volkes forfinancial support and acknowledges support by thestudy program “Biological Physics” of the Elite Net-work of Bavaria. Furthermore, we thank the labora-tory of professor Alexander Bershadsky at WeizmannInsitute of Science in Isreal for providing the REF52cells stably expressing paxillin-YFP.10
EFERENCES REFERENCES
References [1] Yigang Shen, Yaxiaer Yalikun, and Yo Tanaka.Recent advances in microfluidic cell sorting sys-tems.
Sensors and Actuators B: Chemical , 282:268–281, March 2019. ISSN 09254005. doi:10.1016/j.snb.2018.11.025.[2] Oliver Otto, Philipp Rosendahl, Alexander Mi-etke, Stefan Golfier, Christoph Herold, DanielKlaue, Salvatore Girardo, Stefano Pagliara, An-drew Ekpenyong, Angela Jacobi, Manja Wobus,Nicole T¨opfner, Ulrich F Keyser, J¨org Mansfeld,Elisabeth Fischer-Friedrich, and Jochen Guck.Real-time deformability cytometry: On-the-flycell mechanical phenotyping.
Nature Methods ,12(3):199–202, March 2015. ISSN 1548-7091,1548-7105. doi: 10.1038/nmeth.3281.[3] Bob Fregin, Fabian Czerwinski, Doreen Bieden-weg, Salvatore Girardo, Stefan Gross, Kon-stanze Aurich, and Oliver Otto. High-throughputsingle-cell rheology in complex samples by dy-namic real-time deformability cytometry.
Na-ture Communications , 10(1):415, December2019. ISSN 2041-1723. doi: 10.1038/s41467-019-08370-3.[4] Philip Kollmannsberger and Ben Fabry. Linearand Nonlinear Rheology of Living Cells.
AnnualReview of Materials Research , 41(1):75–97, Au-gust 2011. ISSN 1531-7331, 1545-4118. doi:10.1146/annurev-matsci-062910-100351.[5] Rafael D Gonzalez-Cruz, Vera C Fonseca, andEric M Darling. Cellular mechanical propertiesreflect the differentiation potential of adipose-derived mesenchymal stem cells.
Proc. Nat.Acad. Sci. (USA) , 109(24):E1523–E1529, 2012.[6] F Huber, J Schnau, S R¨onicke, P Rauch,K M¨uller, C F¨utterer, and J K¨as. Emergent com-plexity of the cytoskeleton: from single filamentsto tissue.
Advances in Physics , 62(1):1–112,February 2013.[7] Tom Bongiorno, Jacob Kazlow, RomanMezencev, Sarah Griffiths, Rene Olivares-Navarrete, John F McDonald, Zvi Schwartz,Barbara D Boyan, Todd C McDevitt, andTodd Sulchek. Mechanical stiffness as animproved single-cell indicator of osteoblastichuman mesenchymal stem cell differentiation.
JBiomechanics , 47(9):2197–2204, June 2014. [8] Elisabeth Fischer-Friedrich, Anthony A. Hyman,Frank J¨ulicher, Daniel J. M¨uller, and Jonne He-lenius. Quantification of surface tension and in-ternal pressure generated by single mitotic cells.
Scientific Reports , 4:6213, August 2014. ISSN2045-2322. doi: 10.1038/srep06213.[9] Janina R Lange, Julian Steinwachs, ThorstenKolb, Lena A Lautscham, Irina Harder, GraemeWhyte, and Ben Fabry. Microconstriction Ar-rays for High-Throughput Quantitative Measure-ments of Cell Mechanical Properties.
Biophys.J. , 109(1):26–34, July 2015.[10] Elisabeth Fischer-Friedrich, Yusuke Toyoda,Cedric J. Cattin, Daniel J. M¨uller, Anthony A.Hyman, and Frank J¨ulicher. Rheology of the Ac-tive Cell Cortex in Mitosis.
Biophysical Journal ,111(3):589–600, August 2016. ISSN 00063495.doi: 10.1016/j.bpj.2016.06.008.[11] Kendra D Nyberg, Kenneth H Hu, Sara H Klein-man, Damir B Khismatullin, Manish J Butte,and Amy C Rowat. Quantitative DeformabilityCytometry: Rapid, Calibrated Measurements ofCell Mechanical Properties.
Biophys. J. , 113(7):1574–1584, October 2017.[12] Janina R Lange, Claus Metzner, SebastianRichter, Werner Schneider, Monika Spermann,Thorsten Kolb, Graeme Whyte, and Ben Fabry.Unbiased High-Precision Cell Mechanical Mea-surements with Microconstrictions.
Biophys. J. ,112(7):1472–1480, April 2017.[13] H Kubitschke, J Schnauss, K D Nnetu, E Warmt,R Stange, and J Kaes. Actin and microtubulenetworks contribute differently to cell responsefor small and large strains.
New J. Phys. , 19(9):093003–13, September 2017.[14] Devina Jaiswal, Norah Cowley, Zichao Bian,Guoan Zheng, Kevin P. Claffey, and KazunoriHoshino. Stiffness analysis of 3D spheroidsusing microtweezers.
PLoS One , 12(11):e0188346, 2017.[15] Yuval Mulla, F C MacKintosh, and Gijsje HKoenderink. Origin of Slow Stress Relaxationin the Cytoskeleton.
Phys. Rev. Lett. , 122(21):218102, May 2019.[16] Jessica Snyder, Ae Rin Son, Qudus Hamid,Chengyang Wang, Yigong Lui, and Wei Sun.11
EFERENCES REFERENCES
Mesenchymal stem cell printing and processregulated cell properties.
Biofabrication , 7(4):044106, December 2015. ISSN 1758-5090. doi:10.1088/1758-5090/7/4/044106.[17] Andreas Blaeser, Daniela Filipa Duarte Cam-pos, Uta Puster, Walter Richtering, Molly M.Stevens, and Horst Fischer. Controlling ShearStress in 3D Bioprinting is a Key Factor to Bal-ance Printing Resolution and Stem Cell Integrity.
Advanced Healthcare Materials , 5(3):326–333,December 2015. ISSN 2192-2640. doi: 10.1002/adhm.201500677.[18] Yu Zhao, Yang Li, Shuangshuang Mao, WeiSun, and Rui Yao. The influence of printingparameters on cell survival rate and printabilityin microextrusion-based 3D cell printing tech-nology.
Biofabrication , 7(4):045002, Novem-ber 2015. ISSN 1758-5090. doi: 10.1088/1758-5090/7/4/045002.[19] Naomi Paxton, Willi Smolan, Thomas B¨ock,Ferry Melchels, J¨urgen Groll, and TomaszJungst. Proposal to assess printability of bioinksfor extrusion-based bioprinting and evaluation ofrheological properties governing bioprintability.
Biofabrication , 9(4):044107, November 2017.ISSN 1758-5090. doi: 10.1088/1758-5090/aa8dd8.[20] Sebastian J. M¨uller, Elham Mirzahossein,Emil N. Iftekhar, Christian B¨acher, StefanSchr¨ufer, Dirk W. Schubert, Ben Fabry, andStephan Gekle. Flow and hydrodynamic shearstress inside a printing needle during biofabrica-tion.
PLOS ONE , 15(7):e0236371, July 2020.ISSN 1932-6203. doi: 10.1371/journal.pone.0236371.[21] Saif Khalil and Wei Sun. Biopolymer depositionfor freeform fabrication of hydrogel tissue con-structs.
Materials Science and Engineering: C ,27(3):469–478, April 2007.[22] Brian A Aguado, Widya Mulyasasmita, JamesSu, Kyle J Lampe, and Sarah C Heilshorn. Im-proving Viability of Stem Cells During SyringeNeedle Flow Through the Design of HydrogelCell Carriers.
Tissue Engineering Part A , 18(7-8):806–815, April 2012.[23] Annalisa Tirella, Federico Vozzi, GiovanniVozzi, and Arti Ahluwalia. PAM2 (Piston As- sisted Microsyringe): A New Rapid Prototyp-ing Technique for Biofabrication of Cell Incor-porated Scaffolds.
Tissue Engineering Part C:Methods , 17(2):229–237, February 2011.[24] Minggan Li, Xiaoyu Tian, Janusz A Kozinski,Xiongbiao Chen, and Dae Kun Hwang. Mod-eling mechanical cell damage in the bioprintingprocess employing a conical needle.
J. Mech.Med. Biol. , 15(05):1550073–15, October 2015.[25] V. V. Lulevich, I. L. Radtchenko, G. B. Sukho-rukov, and O. I. Vinogradova. DeformationProperties of Nonadhesive Polyelectrolyte Mi-crocapsules Studied with the Atomic Force Mi-croscope.
The Journal of Physical Chemistry B ,107(12):2735–2740, March 2003. ISSN 1520-6106, 1520-5207. doi: 10.1021/jp026927y.[26] Valentin Lulevich, Tiffany Zink, Huan-YuanChen, Fu-Tong Liu, and Gang-yu Liu. Cell Me-chanics Using Atomic Force Microscopy-BasedSingle-Cell Compression.
Langmuir , 22(19):8151–8155, September 2006. ISSN 0743-7463,1520-5827. doi: 10.1021/la060561p.[27] Hamid Ladjal, Jean-Luc Hanus, Anand Pil-larisetti, Carol Keefer, Antoine Ferreira, and Jay-dev P. Desai. Atomic force microscopy-basedsingle-cell indentation: Experimentation and fi-nite element simulation. In , pages 1326–1332, St. Louis, MO, USA,October 2009. IEEE. ISBN 978-1-4244-3803-7.doi: 10.1109/IROS.2009.5354351.[28] Robert Kiss. Elasticity of Human EmbryonicStem Cells as Determined by Atomic Force Mi-croscopy.
Journal of Biomechanical Engineer-ing , 133(10):101009, November 2011. ISSN0148-0731. doi: 10.1115/1.4005286.[29] Fabian M Hecht, Johannes Rheinlaender, Nico-las Schierbaum, Wolfgang H Goldmann, BenFabry, and Tilman E Sch¨affer. Imaging vis-coelastic properties of live cells by AFM: power-law rheology on the nanoscale.
Soft Matter , 11(23):4584–4591, 2015.[30] Ali Ghaemi, Alexandra Philipp, Andreas Bauer,Klaus Last, Andreas Fery, and Stephan Gekle.Mechanical behaviour of micro-capsules andtheir rupture under compression.
Chem. Eng.Sci. , 142(C):236–243, March 2016.12
EFERENCES REFERENCES [31] Ana Sancho, Ine Vandersmissen, Sander Craps,Aernout Luttun, and J¨urgen Groll. A new strat-egy to measure intercellular adhesion forces inmature cell-cell contacts.
Sci. Rep. , 7(1):46152–14, April 2017.[32] Yuri M Efremov, Wen-Horng Wang, Shana DHardy, Robert L Geahlen, and Arvind Raman.Measuring nanoscale viscoelastic parameters ofcells directly from AFM force-displacementcurves.
Sci. Rep. , 7(1):1541–14, May 2017.[33] Hamid Ladjal, Jean-Luc Hanus, Anand Pil-larisetti, Carol Keefer, Antoine Ferreira, and Jay-dev P Desai. Atomic force microscopy-basedsingle-cell indentation: Experimentation and fi-nite element simulation. In , pages 1326–1332. IEEE,September 2018.[34] Ya Hua Chim, Louise M Mason, Nicola Rath,Michael F Olson, Manlio Tassieri, and HuabingYin. A one-step procedure to probe the vis-coelastic properties of cells by Atomic Force Mi-croscopy.
Sci. Rep. , 8(1):1–12, September 2018.[35] Kenneth L. Johnson.
Contact Mechanics . Cam-bridge Univ. Press, Cambridge, 9. print edition,2003. ISBN 978-0-521-34796-9.[36] Edward Dintwa, Engelbert Tijskens, and Her-man Ramon. On the accuracy of the Hertzmodel to describe the normal contact of soft elas-tic spheres.
Granular Matter , 10(3):209–221,March 2008. ISSN 1434-5021, 1434-7636. doi:10.1007/s10035-007-0078-7.[37] Jens W. Neubauer, Nicolas Hauck, Max J.M¨annel, Maximilian Seuss, Andreas Fery, andJulian Thiele. Mechanoresponsive Hydrogel Par-ticles as a Platform for Three-Dimensional ForceSensing.
ACS Applied Materials & Interfaces , 11(29):26307–26313, July 2019. ISSN 1944-8244,1944-8252. doi: 10.1021/acsami.9b04312.[38] Jonathan B Freund. Numerical Simulation ofFlowing Blood Cells.
Annu. Rev. Fluid Mech. ,46(1):67–95, January 2014.[39] G´abor Z´avodszky, Britt van Rooij, Victor Azizi,and Alfons Hoekstra. Cellular Level In-silicoModeling of Blood Rheology with An ImprovedMaterial Model for Red Blood Cells.
Front.Physiol. , 8:061006–14, August 2017. [40] Johannes Mauer, Simon Mendez, Luca Lan-otte, Franck Nicoud, Manouk Abkarian, Ger-hard Gompper, and Dmitry A Fedosov. Flow-Induced Transitions of Red Blood Cell Shapesunder Shear.
Phys. Rev. Lett. , 121(11):118103,September 2018.[41] Achim Guckenberger, Alexander Kihm, ThomasJohn, Christian Wagner, and Stephan Gekle.Numerical-experimental observation of shapebistability of red blood cells flowing in a mi-crochannel.
Soft Matter , 14(11):2032–2043,March 2018.[42] Christos Kotsalos, Jonas Latt, and BastienChopard. Bridging the computational gap be-tween mesoscopic and continuum modeling ofred blood cells for fully resolved blood flow.
J.Comput. Phys. , 398:108905, December 2019.[43] Nathalie Caille, Olivier Thoumine, Yanik Tardy,and Jean-Jacques Meister. Contribution of thenucleus to the mechanical properties of endothe-lial cells.
Journal of Biomechanics , 35(2):177–187, February 2002. ISSN 00219290. doi:10.1016/S0021-9290(01)00201-9.[44] M Mokbel, D Mokbel, A Mietke, N Tr¨aber, S Gi-rardo, O Otto, J Guck, and S Aland. NumericalSimulation of Real-Time Deformability Cytom-etry To Extract Cell Mechanical Properties.
ACSBiomater. Sci. Eng. , 3(11):2962–2973, January2017.[45] R Roscoe. On the rheology of a suspension ofviscoelastic spheres in a viscous liquid.
J. FluidMech. , 28(02):273–21, March 1967.[46] Tong Gao and Howard H. Hu. Deformation ofelastic particles in viscous shear flow.
Journalof Computational Physics , 228(6):2132–2151,April 2009. ISSN 00219991. doi: 10.1016/j.jcp.2008.11.029.[47] Tong Gao, Howard H Hu, and Pedro PonteCasta˜neda. Rheology of a suspension of elasticparticles in a viscous shear flow.
J. Fluid Mech. ,687:209–237, October 2011.[48] Tong Gao, Howard H Hu, and Pedro PonteCasta˜neda. Shape Dynamics and Rheology ofSoft Elastic Particles in a Shear Flow.
Phys. Rev.Lett. , 108(5):058302–4, January 2012.13
EFERENCES REFERENCES [49] Kirill Lykov, Yasaman Nematbakhsh, MenglinShang, Chwee Teck Lim, and Igor V Pivkin.Probing eukaryotic cell mechanics via meso-scopic simulations.
PLoS Comput Biol , 13(9):e1005726–22, September 2017.[50] M M Villone, M A Hulsen, P D Anderson, andP L Maffettone. Simulations of deformable sys-tems in fluids under shear flow using an arbitraryLagrangian Eulerian technique.
Computers &Fluids , 90(C):88–100, February 2014.[51] M M Villone, G D’Avino, M A Hulsen, and P LMaffettone. Dynamics of prolate spheroidal elas-tic particles in confined shear flow.
Phys. Rev. E ,92(6):062303–12, December 2015.[52] Marco E. Rosti, Luca Brandt, and DhrubadityaMitra. Rheology of suspensions of viscoelas-tic spheres: Deformability as an effective vol-ume fraction.
Physical Review Fluids , 3(1):012301, January 2018. ISSN 2469-990X. doi:10.1103/PhysRevFluids.3.012301.[53] Amir Saadat, Christopher J. Guido, GianlucaIaccarino, and Eric S. G. Shaqfeh. Immersed-finite-element method for deformable particlesuspensions in viscous and viscoelastic media.
Physical Review E , 98(6):063316, December2018. ISSN 2470-0045, 2470-0053. doi: 10.1103/PhysRevE.98.063316.[54] Antonina Y. Alexandrova, Katya Arnold,S´ebastien Schaub, Jury M. Vasiliev, Jean-Jacques Meister, Alexander D. Bershadsky,and Alexander B. Verkhovsky. ComparativeDynamics of Retrograde Actin Flow and FocalAdhesions: Formation of Nascent AdhesionsTriggers Transition from Fast to Slow Flow.
PLoS ONE , 3(9):e3234, September 2008. ISSN1932-6203. doi: 10.1371/journal.pone.0003234.[55] Allan F. Bower.
Applied Mechanics of Solids .CRC Press, Boca Raton, 2010. ISBN 978-1-4398-0247-2.[56] M. Mooney. A Theory of Large Elastic Deforma-tion.
Journal of Applied Physics , 11(9):582–592,September 1940. ISSN 0021-8979, 1089-7550.doi: 10.1063/1.1712836.[57] R. S. Rivlin. Large Elastic Deformations ofIsotropic Materials. I. Fundamental Concepts.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sci-ences , 240(822):459–490, January 1948. ISSN1364-503X, 1471-2962. doi: 10.1098/rsta.1948.0002.[58] Christian B¨acher and Stephan Gekle. Computa-tional modeling of active deformable membranesembedded in three-dimensional flows.
Phys. Rev.E , 99(6):062418, June 2019.[59] S. Ramanujan and C. Pozrikidis. Deforma-tion of liquid capsules enclosed by elastic mem-branes in simple shear flow: Large deforma-tions and the effect of fluid viscosities.
Journalof Fluid Mechanics , 361:117–143, April 1998.ISSN 0022-1120, 1469-7645. doi: 10.1017/S0022112098008714.[60] Jonathan R. Clausen and Cyrus K. Aidun. Cap-sule dynamics and rheology in shear flow: Parti-cle pressure and normal stress.
Physics of Fluids ,22(12):123302, December 2010. ISSN 1070-6631, 1089-7666. doi: 10.1063/1.3483207.[61] Achim Guckenberger, Marcel P Schraml, Paul GChen, Marc Leonetti, and Stephan Gekle. Onthe bending algorithms for soft objects in flows.
Comput. Phys. Commun. , 207:1–23, October2016.[62] Michael Smith.
ABAQUS/Standard User’s Man-ual, Version 6.9 . Dassault Syst`emes SimuliaCorp, United States, 2009.[63] H.J. Limbach, A. Arnold, B.A. Mann, andC. Holm. ESPResSo—an extensible simula-tion package for research on soft matter sys-tems.
Computer Physics Communications , 174(9):704–727, May 2006. ISSN 00104655. doi:10.1016/j.cpc.2005.10.005.[64] D. Roehm and A. Arnold. Lattice Boltzmannsimulations on GPUs with ESPResSo.
The Eu-ropean Physical Journal Special Topics , 210(1):89–100, August 2012. ISSN 1951-6355, 1951-6401. doi: 10.1140/epjst/e2012-01639-6.[65] Dharshi Devendran and Charles S Peskin. Animmersed boundary energy-based method for in-compressible viscoelasticity.
J. Comput. Phys. ,231(14):4613–4642, May 2012.[66] Christian B¨acher, Lukas Schrack, and StephanGekle. Clustering of microscopic particles in14
EFERENCES REFERENCES constricted blood flow.
Phys. Rev. Fluids , 2(1):013102, January 2017.[67] Moritz Lehmann, Sebastian Johannes M¨uller,and Stephan Gekle. Efficient viscosity contrastcalculation for blood flow simulations using thelattice Boltzmann method.
Int. J. Numer. Meth.Fluids , 103(18):1–15, April 2020. 15 S UPPLEMENTARY MATERIAL FOR THE MANUSCRIPT
A hyperelastic model for simulating cells in flow
Sebastian J. M¨uller , Franziska Weigl , Carina Bezold , Ana Sancho , , ChristianB¨acher , Krystyna Albrecht , and Stephan Gekle Theoretical Physics VI, Biofluid Simulation and Modeling, University ofBayreuth, Universit¨atsstraße 30, 95440 Bayreuth, Germany Department of Functional Materials in Medicine and Dentistry and BavarianPolymer Institute (BPI), University of W¨urzburg, Pleicherwall 2, 97070 W¨urzburg,Germany Department of Automatic Control and Systems Engineering, University of theBasque Country UPV/EHU, San Sebastian, SpainAuthor to whom correspondence should be addressed: [email protected]
S-1 Supplementary Material for the cell experiments
Additional force-deformation curves for our FluidFM R (cid:13) measurements on REF52cells are shown in figure S-1. Compared to the curves depicted in the manuscript infigure 4, these measurements show an earlier upturn of the force. Thus, our modeloverestimates the force necessary for a small deformation of the cell and slightly un-derestimates the force for larger deformations. Nevertheless, all measurements fit inthe simulated range of E = ±
100 Pa for w = .
25 and an averaged cell radius of8 . ( ) µm, as figure S-1 shows. The cell radii and Young’s moduli for all measure-ments are listed in table S-1. Table S-1
Measured cell radii R and fitted Young’s moduli E and w for our FluidFM R (cid:13) experiments.Number 1 2 3 4 5 6 7 8 9 R [ µm ] E [ Pa ]
160 190 220 170 210 290 210 220 125 w N o r m a l f o r ce [ n N ] Deformation [%]
Experiment Our model56789 210 Pa290 Pa210 Pa220 Pa125 Pa
Fig. S-1
Our numerical model in comparison to our FluidFM R (cid:13) measurements on REF52 cells. The ratioof the shear moduli is chosen as w = .
25 for all curves. The gray area shows the simulation of a cell withan averaged cell radius of 8 . ( ) µm and Young’s modulus range 220 ± S-2 Supporting Information for the numerical model
S-2.1 Convergence of single cell deformation in shear flowThe temporal development of the deformation D of a single cell in a Couette flowcan be seen in figure S-2. Starting from a spherical shape ( D = > .
2, we first find an overrelaxation of thedeformation before it converges towards a constant value. . . . . . . . . . D step / − Ca = . = . = . = . = . Fig. S-2
Single cell deformation in Couette flow for different capillary numbers. After an initial transienttimespan, the deformation converges to a constant value.
S-2.2 Reduction of the system resolutionIn figure S-3 we show that a system with reduced cell resolution (from R Cell =
10 to R Cell = × ×
100 to 60 × ×
30 grid cells) produces the same deformation versus capillary number behavior as thesystem with higher resolution.S-2.3 Translational and rotational invariance of the force calculationAs a very direct test for the correct behavior of our model, we consider a single tetra-hedron and examine the behavior of the volume and the elastic force for an initiallyapplied translation, rotation and stretching. In figure S-4a, the behavior of the vol-ume under these deformations is shown over the first time steps. While the volumeremains constant under pure translation, pure rotation, and a combination of both,it quickly relaxes towards its reference value after an initial stretch is applied. Thesame behavior is observed for the elastic force acting on one tetrahedron vertex, infigure S-4b. . . . . . . . . D Calow resolutionhigh resolution
Fig. S-3
Taylor deformation as function of the capillary number for two different cell and channel reso-lutions. The large system ( R Cell =
10, box: 100 × ×
100 grid cells) produces the same outcome as thedown-scaled system ( R Cell =
6, box: 60 × ×
30 grid cells).(a) . . . . . .
06 0 20 40 60 80 100 V o l u m e / r e f e r e n ce v o l u m e Integration step pure rotationpure translationrotation + translationpure stretchstretch + rotationstretch + translationstretch + rotation + translation (b) . . . . . . .
14 0 20 40 60 80 100 F o r ce ( m ag n i t ud e ) Integration step pure rotationpure translationrotation + translationpure stretchstretch + rotationstretch + translationstretch + rotation + translation
Fig. S-4
The behavior of (a) the volume and (b) the elastic force on a single vertex of a tetrahedron afteran initial rotation, translation or stretching.
S-2.4 Mesh generation and mesh independenceThe tetrahedral mesh of our spheroid is generated using the software gmsh (version4 . .
0) [1]. The Frontal2D meshing algorithm produced a mesh with highest unifor-mity considering edge length, triangle area and tetrahedron volume distribution. Nev-ertheless, all other available meshing algorithms produce likewise uniform meshes,with one exception being the Frontal3D algorithm, as listed in table S-2. We de-mand the uniformity of the mesh to increase the accuracy of our coupled Immersed-Boundary Lattice Boltzmann simulations. Figure S-5 shows the force-deformationcurves for meshes with increasing number of tetrahedra, which are converged andthus prove sufficient sampling of the volume mesh.
Table S-2
Statistics of meshes created using different built-in algorithms of Gmsh [1]. Listed are edgelength L , triangle area A , and tetrahedron volume V providing average, standard deviation, minimum andmaximum value for each mesh.Algorithm Frontal2D MeshAdapt Delaunay2D Delaunay3D Frontal3D¯ L σ L L min L max A σ A A min A max V σ V V min V max N o r m a l f o r ce [ n N ] Deformation [%]460733586624688 ? † ‡ ? used for simulations † sufficiently converged ‡ sufficiently stable Fig. S-5
Force-deformation behavior of meshes with increasing number of tetrahedra. Meshes with N ≥ R = . E = ν = . S-2.5 Hertz theoryAlthough originally designed for the contact between two linear elastic spheres, theHertz theory can be applied to the contact between a linear elastic sphere and a flatplate [2]. The general assumptions for the Hertz-theory are the following [3, p. 91-92]: – frictionless, smooth contact surfaces – contact area small compared to sphere dimension – homogeneous, isotropic and linear elastic material S-2.5.1 Sphere-sphere contact
The following quantities are necessary to describe the normal contact of two elasticspheres. The radii R and R of the spheres define the effective radius of curvature R of the bodies by 1 R = R + R . (S-1)Through their Young’s moduli and the Poisson ratios, E , E and ν , ν , the effectivestiffness K is defined as: 1 K = − ν E + − ν E (S-2)The displacement δ , which measures the distance that the sphere centers approacheach other due to a normal force N acting on each sphere, can be expressed in termsof the above parameters [2]: δ = (cid:18) N RK (cid:19) (S-3)Therefore, the force–displacement relation according to the Hertzian theory for asphere-sphere contact is given by N ( δ ) = KR δ . (S-4) S-2.5.2 Sphere-plane contact
The analytical solution for the force–displacement relation according to the Hertziantheory for the contact of a linear elastic sphere with a rigid plane can be obtainedfrom (S-4) by applying the following modifications: the plane has no curvature, thus R → ∞ and (S-1) simply yields R = R . Since the plane is assumed rigid, i. e. E (cid:29) E , (S-2) reduces to K = E − ν . In this case, N is the force acting on the sphere and δ is the distance between the center of the sphere and the plane.S-2.6 Influence of the Mooney-Rivlin ratio w To clarify the influence of w , we plot in figure S-6 the force versus deformation be-havior of our cell model for different values of w . With decreasing w , i. e. decreasing µ while increasing µ , the strain hardening effect clearly increases and the upturn ofthe force curve begins at lower deformations. This is due to µ scaling the term in thestrain energy density that is quadratic with the deformation (cf. equations (4) and (5)of the manuscript).S-2.7 Influence of the Poisson ratio ν In figure S-7 we demonstrate that variations of the Poisson ratio ν within the rangeof an approximately incompressible material do not notably influence the force-deformation curves. N o r m a l f o r ce [ n N ] Deformation [%] w = 1 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . Fig. S-6
Variation of w : The lower w (the lower µ compared to µ ), the stronger the non-linear upturn ofthe force becomes. The curve with w = E = N o r m a l f o r ce [ n N ] Deformation [%] ν = 0 . ν = 0 . ν = 0 . ν = 0 . ν = 0 . ν = 0 . Fig. S-7
Force versus deformation curves for different Poisson ratios ν . The following parameters wereused: cell radius R = . E = S-3 Compression and indentation simulations
After initialization, each time step of our overdamped relaxation simulation consistsof the following two steps: the movement of the upper wall to compress – or thesphere to indent – the cell and the integration of the equation of motion of the cellvertices, ˙ y α = γ − ( f α + f α probe ) . (S-5)The vertex velocity ˙ y α is obtained from the elastic restoring forces ( f α (12) and theprobe repulsion f α probe ) , considering a friction factor γ . Since here we are only look-ing at a sequence of equilibrium states, the value of γ is irrelevant for the resultingforce-deformation curves and only influences the performance and stability of thesimulations. The equation of motion is integrated using a fourth order Runge-Kuttaalgorithm. The repulsive cell-probe interaction, preventing the cell vertices from pen-etrating the plates or the indenter, has the form f probe ( d ) = c F d n , (S-6)with the cell-probe distance d and a proportionality factor c F . The force points normalto the probe, resulting in a compression between two plates and a radial displacementaway from the indenter. Physically, this corresponds to a free-slip boundary conditionwhich does not restrict tangential motions of the cell along the probe. S-4 Flow simulations with Lattice Boltzmann
S-4.1 MethodThis section briefly summarizes the Lattice Boltzmann method implemented in theopen-source package ESPResSo [4]. For an introduction into the Lattice Boltzmannmethod we refer the interested reader to the book by Kr¨uger et al. [5]. The LatticeBoltzmann equation for the multiple relaxation time scheme used in ESPResSo reads: f i ( x + c i ∆ t , t + ∆ t ) − f i ( x , t ) = ∑ j = (cid:0) M − ω M (cid:1) i j (cid:16) f j ( x , t ) − f eq j ( x , t ) (cid:17) (S-7)It describes the collision and streaming of the population distribution f i ( i = , . . . , ∆ t . Here, c i are the discretized lattice velocities, M denotes trans-formation matrix that maps the populations onto moment space, ω is the diagonalrelaxation frequency matrix, and f eq i denote the equilibrium population distributions.The relaxation frequency for the shear moments ω S is related to the dynamic viscosityof the fluid via [6] η = ρ c s (cid:18) ω S − (cid:19) ∆ t , (S-8)with the fluid mass density ρ and the lattice speed of sound c s . In order to ensuresimulation stability, we choose the time step globally according to Kr¨uger et al. [5,p. 273] as ∆ t = c s (cid:0) τ − (cid:1) ∆ x ν ˜ t = ∆ x ν ˜ t , (S-9)with c s = , a global relaxation parameter τ =
1, the kinematic viscosity ν , and anadditional factor ˜ t in the range 1–2 to manually tune the time step.We further introduce a scaling factor r by which we divide both the viscosity and theYoung’s modulus. According to eq. (S-9), this leads to a larger time step and thus to aspeed-up of the simulations. At the same time it leaves the important Capillary num-ber unchanged and only increases the Reynolds number, which nevertheless remains (cid:28)
1. The parameter r thus does not affect the physics of the simulation which wehave carefully checked by a number of test runs with r = × ×
30 ( x × y × z ), for themultiple cell simulation (cf. section 6.2) it is 50 × ×
40. The dynamic viscosity,chosen as ν = ∆ t = with ˜ t = − − − − − − − − − − − y / R Ca = 0 . . y / R x/R Ca = 1 . x/R Ca = 1 . Fig. S-8
The trajectory of a surface node (here: starting at y = R and x , z =
0) for different capillarynumbers traces the ellipsoidal contour of the deformed particle. The non-elliptical part of the trajectory inthe upper-right corner represents the approach from the initially spherical to the final shape.
S-4.2 Tank-treading motionFigure S-8 shows the trajectories of selected vertices on the outer surface of the par-ticle for different capillary numbers. They describe an ellipsoidal motion tracing theouter contour of the deformed particle thus demonstrating that in our simulations theparticle exhibits tank-treading.
References
1. C. Geuzaine, J.F. Remacle, International Journal for Numerical Methods in Engi-neering (11), 1309 (2009). DOI 10.1002/nme.25792. E. Dintwa, E. Tijskens, H. Ramon, Granular Matter (3), 209 (2008). DOI10.1007/s10035-007-0078-73. K.L. Johnson, Contact Mechanics , 9th edn. (Cambridge Univ. Press, Cambridge,2003). OCLC: 2500043674. H. Limbach, A. Arnold, B. Mann, C. Holm, Computer Physics Communications (9), 704 (2006). DOI 10.1016/j.cpc.2005.10.0055. T. Kr¨uger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen,
TheLattice Boltzmann Method . Graduate Texts in Physics (Springer International Pub-lishing, Cham, 2017). DOI 10.1007/978-3-319-44649-3
6. Z. Chai, B. Shi, Z. Guo, F. Rong, Journal of Non-Newtonian Fluid Mechanics166