A hyperelastic model for simulating cells in flow

In the emerging field of 3D bioprinting, cell damage due to large deformations is considered a main cause for cell death and loss of functionality inside the printed construct. Those deformations, in turn, strongly depend on the mechano-elastic response of the cell to the hydrodynamic stresses experienced during printing. In this work, we present 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 employ force calculations on a tetrahedralized volume mesh. To calibrate our model, we perform a series of FluidFM\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{{\textregistered }}$$\end{document}® compression experiments with REF52 cells demonstrating that all three parameters of the Mooney–Rivlin model are required for a good description of the experimental data at very large deformations up to 80%. In addition, we validate the model by comparing to previous AFM experiments on bovine endothelial cells and artificial hydrogel particles. To investigate cell deformation in flow, we incorporate our model into Lattice Boltzmann simulations via an Immersed-Boundary algorithm. In linear shear flows, our model shows excellent agreement with analytical calculations and previous simulation data. Electronic supplementary material The online version of this article (10.1007/s10237-020-01397-2) contains supplementary material, which is available to authorized users.


S-1 Supplementary Material for the cell experiments
Additional force-deformation curves for our FluidFM R measurements on REF52 cells are shown in figure S-1. Compared to the curves depicted in the manuscript in figure 4, these measurements show an earlier upturn of the force. Thus, our model overestimates the force necessary for a small deformation of the cell and slightly underestimates the force for larger deformations. Nevertheless, all measurements fit in the simulated range of E = 220 ± 100 Pa for w = 0.25 and an averaged cell radius of 8.6(7) µm, as figure S-1 shows. The cell radii and Young's moduli for all measurements are listed in table S-1.  The gray area shows the simulation of a cell with an averaged cell radius of 8.6(7) µm and Young's modulus range 220 ± 100 Pa.

S-2.1 Convergence of single cell deformation in shear flow
The temporal development of the deformation D of a single cell in a Couette flow can be seen in figure S-2. Starting from a spherical shape (D = 0), the cell experiences a shape change during an initial transient timespan, after which it assumes a steady shape. For capillary numbers Ca > 0.2, we first find an overrelaxation of the deformation before it converges towards a constant value.  S-2.2 Reduction of the system resolution In figure S-3 we show that a system with reduced cell resolution (from R Cell = 10 to R Cell = 6 grid cells) and a smaller simulation box (from 100 × 150 × 100 to 60 × 90 × 30 grid cells) produces the same deformation versus capillary number behavior as the system with higher resolution.

S-2.3 Translational and rotational invariance of the force calculation
As a very direct test for the correct behavior of our model, we consider a single tetrahedron and examine the behavior of the volume and the elastic force for an initially applied translation, rotation and stretching. In figure S-4a, the behavior of the volume under these deformations is shown over the first time steps. While the volume remains constant under pure translation, pure rotation, and a combination of both, it quickly relaxes towards its reference value after an initial stretch is applied. The same behavior is observed for the elastic force acting on one tetrahedron vertex, in figure S-4b.

S-2.4 Mesh generation and mesh independence
The tetrahedral mesh of our spheroid is generated using the software gmsh (version 4.3.0) [1]. The Frontal2D meshing algorithm produced a mesh with highest uniformity considering edge length, triangle area and tetrahedron volume distribution. Nevertheless, all other available meshing algorithms produce likewise uniform meshes, with one exception being the Frontal3D algorithm, as listed in table S-2. We demand the uniformity of the mesh to increase the accuracy of our coupled Immersed-Boundary Lattice Boltzmann simulations. Figure S-5 shows the force-deformation curves for meshes with increasing number of tetrahedra, which are converged and thus prove sufficient sampling of the volume mesh.

S-2.5 Hertz theory
Although originally designed for the contact between two linear elastic spheres, the Hertz theory can be applied to the contact between a linear elastic sphere and a flat plate [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 elastic spheres. The radii R 1 and R 2 of the spheres define the effective radius of curvature R of the bodies by Through their Young's moduli and the Poisson ratios, E 1 , E 2 and ν 1 , ν 2 , the effective stiffness K is defined as: The displacement δ , which measures the distance that the sphere centers approach each other due to a normal force N acting on each sphere, can be expressed in terms of the above parameters [2]: Therefore, the force-displacement relation according to the Hertzian theory for a sphere-sphere contact is given by .

S-2.5.2 Sphere-plane contact
The analytical solution for the force-displacement relation according to the Hertzian theory for the contact of a linear elastic sphere with a rigid plane can be obtained from (S-4) by applying the following modifications: the plane has no curvature, thus R 2 → ∞ and (S-1) simply yields R = R 1 . Since the plane is assumed rigid, i. e. E 2 E 1 , (S-2) reduces to K = E 1 1−ν 2 1 . 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 behavior of our cell model for different values of w. With decreasing w, i. e. decreasing µ 1 while increasing µ 2 , the strain hardening effect clearly increases and the upturn of the force curve begins at lower deformations. This is due to µ 2 scaling the term in the strain 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 range of an approximately incompressible material do not notably influence the forcedeformation curves.

S-3 Compression and indentation simulations
After initialization, each time step of our overdamped relaxation simulation consists of the following two steps: the movement of the upper wall to compress -or the sphere to indent -the cell and the integration of the equation of motion of the cell vertices,ẏ The vertex velocityẏ α is obtained from the elastic restoring forces (f α (12) and the probe repulsion f α probe ), considering a friction factor γ. Since here we are only looking at a sequence of equilibrium states, the value of γ is irrelevant for the resulting force-deformation curves and only influences the performance and stability of the simulations. The equation of motion is integrated using a fourth order Runge-Kutta algorithm. The repulsive cell-probe interaction, preventing the cell vertices from penetrating the plates or the indenter, has the form with the cell-probe distance d and a proportionality factor c F . The force points normal to the probe, resulting in a compression between two plates and a radial displacement away from the indenter. Physically, this corresponds to a free-slip boundary condition which does not restrict tangential motions of the cell along the probe.

S-4.1 Method
This section briefly summarizes the Lattice Boltzmann method implemented in the open-source package ESPResSo [4]. For an introduction into the Lattice Boltzmann method we refer the interested reader to the book by Krüger et al. [5]. The Lattice Boltzmann equation for the multiple relaxation time scheme used in ESPResSo reads: It describes the collision and streaming of the population distribution f i (i = 0, . . . , 18) during one time step ∆t. Here, c i are the discretized lattice velocities, M denotes transformation matrix that maps the populations onto moment space, ω is the diagonal relaxation frequency matrix, and f eq i denote the equilibrium population distributions. The relaxation frequency for the shear moments ω S is related to the dynamic viscosity of the fluid via [6] η = ρc 2 with the fluid mass density ρ and the lattice speed of sound c s . In order to ensure simulation stability, we choose the time step globally according to Krüger et al. [5, p. 273] as with c 2 s = 1 3 , a global relaxation parameter τ = 1, the kinematic viscosity ν, and an additional factort 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 the Young's modulus. According to eq. (S-9), this leads to a larger time step and thus to a speed-up of the simulations. At the same time it leaves the important Capillary number unchanged and only increases the Reynolds number, which nevertheless remains 1. The parameter r thus does not affect the physics of the simulation which we have carefully checked by a number of test runs with r = 1. At the boundaries of the channel a bounce-back algorithm is applied to realize a no-slip boundary condition. For the plane Couette setup, the bounce-back algorithm additionally allows for a fixed tangential velocity component. We use a combined CPU/GPU implementation which enables the calculation of the flow field on the GPU, while the calculation of the cell motion is done in parallel on multiple (4 to 20) CPUs. In lattice units, our simulation box for the single cell in shear flow setup (cf. section 6.1) has the dimensions 60 × 90 × 30 (x × y × z), for the multiple cell simulation (cf. section 6.2) it is 50 × 80 × 40. The dynamic viscosity, chosen as ν = 1 in simulation units, determines the time step in our simulations as ∆t = 1 3 witht = 2. They describe an ellipsoidal motion tracing the outer contour of the deformed particle thus demonstrating that in our simulations the particle exhibits tank-treading.