Numerical Simulation of Real-Time Deformability Cytometry to Extract

Subscriber access provided by UB + Fachbibliothek Chemie | (FU-Bibliothekssystem)

Article

Numerical Simulation of Real-Time Deformability Cytometry to Extract Cell Mechanical Properties Marcel Mokbel, Dominic Mokbel, Alexander Mietke, Nicole Träber, Salvatore Girardo, Oliver Otto, Jochen Guck, and Sebastian Aland ACS Biomater. Sci. Eng., Just Accepted Manuscript • DOI: 10.1021/acsbiomaterials.6b00558 • Publication Date (Web): 11 Jan 2017 Downloaded from http://pubs.acs.org on January 19, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

ACS Biomaterials Science & Engineering is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Biomaterials Science & Engineering

Numerical Simulation of Real-Time Deformability Cytometry to Extract Cell Mechanical Properties M. Mokbel,† D. Mokbel,†,‡ A. Mietke,¶,§ N. Träber,‡ G. Salvatore,‡ O. Otto,‡,k J. Guck,‡ and S. Aland∗,†,⊥ †Institute of Scientific Computing, TU Dresden, Zellescher Weg 12-14, 01069 Dresden, Germany ‡Biotechnology Center, TU Dresden, Tatzberg 47-49, 01307 Dresden, Germany ¶Max-Planck-Institute for Cell Biology and Genetics, Pfotenhauerstrasse 108, 01307 Dresden, Germany §Max-Planck-Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany kCenter for Innovation Competence, University of Greifswald, Fleischmannstrasse 42-44, 17489 Greifswald, Germany ⊥Faculty of Informatics/Mathematics, HTW Dresden, Friedrich-List-Platz 1, 01069 Dresden, Germany E-mail: [email protected]

Abstract The measurement of cell stiffness is an important part of biological research with diverse applications in biology, biotechnology and medicine. Real-time deformability cytometry (RTDC) is a new method to probe cell stiffness at high throughput by flushing cells through a microfluidic channel where cell deformation provides an indicator for cell stiffness. 1 Here,

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

we propose a full numerical model for single cells in a flow channel to quantitatively relate cell deformation to mechanical parameters. Thereby the cell is modeled as a visco-elastic material surrounded by a thin shell cortex, subject to bending stiffness and cortical surface tension. For small deformations our results show good agreement with a previously developed analytical model that neglects the influence of cell deformation on the fluid flow. 2 Including linear elasticity as well as neo-Hookean hyperelasticity, our model is valid in a wide range of cell deformations and allows to extract cell stiffness for largely deformed cells. We introduce a new measure for cell deformation that is capable to distinguish between deformation effects stemming from cell cortex and cell bulk elasticity. Finally, we demonstrate the potential of the method to simultaneously quantify multiple mechanical cell parameters by RT-DC.

Keywords Finite-Element simulation, RT-DC, cell mechanics, elastic moduli, cell stiffness

1

Introduction

The stiffness of cells is a powerful biophysical marker for cell state. Information on cell elasticity can, for instance, be used to distinguish between different cell phenotypes, or between healthy and diseased cells. 3–7 Therefore, the measurement of a cell’s elastic response is an important part of biological research including medical applications in cell sorting and diagnostics 8–16 The heterogeneity of cells in medical and biological samples often enables meaningful statistical analysis of cell mechanical properties only for larger population of cells. Recently, several microfluidic techniques have been introduced that achieve the required high-throughput measurement rates. 1,16–18 One of these methods is real-time deformability cytometry (RT-DC), 1 where suspended animal cells are advected by a shear flow through a microfluidic channel at a constant speed. In this process, cells are deformed due to strong velocity gradients within the channel

2


Page 2 of 34

Page 3 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cross section. 1 Shortly after channel entry, cells assume a stationary shape which is imaged by a high-speed camera. A sketch of the measurement setup is depicted in Fig. 1 (left). The observed deformation of initially spherical cells in the flow channel depends on cell stiffness, but is also influenced by flow speed and relative cell size. Additionally, the internal cell architecture plays an important role and it is currently unknown to which extent the cell bulk and the cortex each contribute to the effective mechanical response of the cell. To disentangle mutual contributions of cell size and cell stiffness to cell deformation an analytic study was presented by Mietke et al. 2 There, the analytical solution for the flow around a sphere within an axisymmetric channel is used to calculate the corresponding surface stresses and to predict the expected shape deformation assuming elastic shell or bulk properties. A comparison between predicted and experimentally measured cell deformations permits conclusions on the elastic moduli of the corresponding cells. Because the analytical model neglects the interaction between the deformed cell shapes and changes in the surrounding hydrodynamic flow profiles, it is only valid for small deformations. Still, the theoretical results of Mietke et al. 2 are regarded as a highly important step toward combining speed with precision in cell mechanics, which can open the door to high-throughput screening and sorting of large cell populations based on the precise mechanical properties of each individual cell. 19 The goal of the present paper is to go beyond the analytical approach in Mietke et al. 2 by introducing a full numerical model that is valid in a wider parameter range. This will allow to probe and analyze the limits of validity of the analytical model in more detail and to extract cell mechanical parameters for stronger deformed cells at finite Reynolds number. To achieve this, the model will include linear elasticity as well as a neo-Hookean hyperelastic model for both, the cell bulk as well as the cell membrane, combined with a cortical tension model. The neo-Hookean model is based on the statistical thermodynamics of cross-linked polymer chains which roughly correspond to the cell cortex and cytoskeleton. The neo-Hookean hyperelastic model is known to be reasonably accurate when the maximum strain is on the order of 100% 20 and has been used for incompressible elasticity of cell membranes, 20,21 cell bulk 22 and whole tissues. 23,24

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Current numerical methods for biological cells in flow are mainly designed for red blood cells (RBCs), which are relatively soft due to the absence of nucleus and most interior organelles. Immersed boundary (IB) methods are widely used for such cells and can naturally include membrane elasticity 21 as well as surface tension and bending stiffness. 25,26 Very recently IB methods have been extended to incompressible visco-elasticity 27 and to coupled membrane and bulk elasticity. 28 Main problems of the IB method are the handling of variable viscosity, 29,30 and the stiffness occurring from the coupling between flow and membrane advection. Efficient parallelization of an adaptive grid IB method is not straightforward and continues to be the subject to active research. 31 An alternative approach are particle methods. Early attempts of dissipative particle dynamics account for cell bulk elasticity and simulate RBCs as solid elastic bodies. 32 Recently, a lot of studies on RBC dynamics used the multiparticle collision dynamics model, where cells are described as elastic capsules with bending rigidity. This numerical scheme is very flexible and allows to simulate many-body problems of vesicles and red blood cells, 33 however the model is partly phenomenological and has to our knowledge not yet been extended to include cell bulk elasticity. Boundary Element Methods (BEM) as used by Pozrikidis 34 are another approach, which can couple the Stokes equations to thin shell theory for the cell membrane. Additional cell bulk elasticity has not been included in the BEM framework. All of the above methods belong to the class of interface tracking methods where an explicit representation of the interface by individual points or a grid is used. In contrast to this, interface capturing methods describe the interface implicitly by an auxiliary field. Elastic contributions are traditionally not included in such methods, although very recently some interesting progress in this direction has been made for bulk elasticity in phase-field, 35 Level-Set 36 and Volume-Of-Fluid methods. 37 Grid based fluid-structure interaction as in the standard benchmark 38 is usually not considered for flowing cells, since its limitation to small displacements. However, in the special case of RT-DC a co-moving grid can be employed to keep the cell in the center of the computational domain. The

4


Page 4 of 34

Page 5 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


representation of the fluid and cell domain by two fitted grids enables a very accurate description of the cell surface which makes the approach very efficient. In contrast to most of the discussed numerical methods, cell bulk elasticity, membrane elasticity and cortical surface tension can all be easily included in one model. We will follow this idea and present an axisymmetric body-fitted numerical approach for single cells in a flow channel. The numerical results will be validated with the analytical results from Mietke et al. 2 for small deformations and can be used to analyze the errors that are made in there due to the neglect of cell deformation on the fluid flow. Using hyperelastic material laws the numerical model can be used in RT-DC to extract elastic parameters for the larger deformations, that are experimentally observed. Comparison to experiments will show that experimental cell shapes can be reproduced in particular if the cell is assumed to be an elastic shell. By combination of bulk and membrane elasticity the model allows to discriminate between both effects and to extract both elastic parameters from cell deformation images.

2

Mechanical Cell Model

The basic structure of eukaryotic cells is to a large extent conserved across different types of animal cells. An impermeable lipid bilayer membrane surrounds an intracellular region filled with fluid and organelles. The membrane is attached to the cell cortex, a thin polymer network at the cell surface that behaves elastically at short timescales and is subject to surface tension forces arising from pulling myosin motor proteins. 39 Coarse graining the elasticity imposed by intracellular organelles and cytoskeleton leads to a representation of the cell as a visco-elastic bulk material, 32 enclosed by a membrane subject to shear and stretch elasticity, bending stiffness and surface tension. Fig. 1 (right) illustrates this picture of a cell and the corresponding forces. We focus on the simulations of single cells on the time scale of milliseconds, in agreement with the time span of a cell in an RT-DC channel. While on larger timescales (e.g., minutes) cells can dissipate elastic stresses by changing their internal structure, this is should not be possible within

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 6 of 34

Page 7 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where vi , Si , ρi are the velocity fields, stresses and densities in both phases, respectively. Note, that the velocity fields are advected relative to the cell velocity, due to the co-moving domains. Viscous and elastic bulk stresses are defined as S0 = −p0 I + η0 ∇v0 + (∇v0 )T , S1 = −p1 I + η1 ∇v1 + (∇v1 )T + Selastic where pi , ηi , I are pressures, viscosities and identity matrix, respectively. The elastic contribution in the cell is described by two different models for incompressible elasticity. The Euler-Almansi strain is the constitutive law underlying the analytical model in Mietke et al. 2 Although it is geometrically non-linear, we will also refer to this model as linear elastic model, since it is the standard model of linear elasticity theory with linear relationships between the components of stress and strain. The second model we employ is the hyperelastic neo-Hookean model which is known to be reasonably accurate when the maximum strain is on the order of 100%. 20 Accordingly, we define the elastic stress,

Selastic =

          

E 3

∇u1 + ∇uT1 − E3 ∇u1 ∇uT1 (linear) E (B 3

− I)

(3)

(neo-Hookean)

where E is the Young’s modulus of the cell bulk, u1 is the displacement field within Ω1 and B is the left Cauchy-Green deformation tensor. The latter two quantities are obtained from the evolution equations,

∂t u1 + (v1 − vcell ) · ∇u1 = v1

in Ω1 .

(4)

∂t B + (v1 − vcell ) · ∇B = (∇v1 )T B + B(∇v1 )

in Ω1 .

(5)

where u1 = 0, B = I prescribe the undeformed reference configuration as initial condition. 35 At

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 34

the interface Γ the following jump conditions are specified to ensure continuity of velocities and balance of forces

[v]Γ = 0,

[S]Γ · n = −

δEtension δEbend δEstretch − − , δΓ δΓ δΓ

where the interface normal n is defined to point into Ω1 and the jump operator is [f ]Γ = f1 − f0 . The interfacial forces are given by the first variation of the interfacial energies with respect to changes in Γ. The surface tension energy and bending energy are given by

Etension =

Z

γ dA,

Ebend =

Γ

Z

Γ

cb 2 κ dA, 8

where cb is the joint bending stiffness of the bilayer and cell cortex, γ is the effective surface tension due to myosin contraction and κ = ∇ · n is the total curvature (twice the mean curvature). The resulting surface tension and bending forces are, 40 δEtension = −γκn, δΓ

δEbend cb = (2∆Γ κ + κ3 − 4κK)n, δΓ 4

where K is the Gaussian curvature of the cell surface. The mathematical definition of the stretching shell energy Estretch is complicated and in general involves tangential surface tensors. 41 However, due to the spherical ground state of cells immersed in a fluid, it is possible to restrict computations to axisymmetric scenarios which leads not only to an enormous reduction in computational complexity but also to a simplified stretching energy. The axisymmetric stretching energy can be expressed in terms of the two principal stretches, that describe the relative change of a surface length, in lateral and rotational direction, respectively, as illustrated in Fig. 10 in the appendix,

λ1 =

ds , ds0

λ2 =

8

r , r0


(6)

Page 9 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where s is the arc length, r the distance to the symmetry axis and s0 , r0 are the corresponding quantities at the same material point in the undeformed state. The cortex stretch elasticity can now be formulated as 42,43

Estretch =

Z

(7)

f (λ1 , λ2 ) dA. Γ

with

f=

          

E2D 1−ν 2

E2D 6

1 (λ1 2

− 1)2 + 21 (λ2 − 1)2 + ν(λ1 − 1)(λ2 − 1)

−2 λ21 + λ22 + λ−2 1 λ2 − 3

(thin shell) ,

(8)

(neo-Hookean shell)

where E2D is the elastic surface modulus and Poisson’s ratio ν = 0.5 for an incompressible shell. The interfacial force resulting from the shell stretching energy is derived in the Appendix. If a square channel is employed in an RT-DC setup, the flow can be mapped onto a cylindrical channel flow using the so-called equivalent channel radius. This radius R is chosen such that there is the same pressure drop along the cylindrical channel and the square channel of width L, which leads to the relation R = 0.547L. 44 It has been shown that this methodology allows accurate axisymmetric predictions of flow speed of spherical objects in square channels as long as the cell occupies less than 90% of the circular channel cross section. 2 The corresponding surface stresses exerted by the fluid are accurately predicted on the lateral sides of the object, but overestimated at the diagonal sides of the object due to the increased distance to the channel corners in a square channel. It is currently not clear to which extent these additional stress inhomogeneities impinge on the cell deformation. However, the good agreement of the circular channel simulations with experimental measurements of elastic beads suggests this is a minor effect for moderate deformations and cell sizes. 2

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3

Page 10 of 34

Discretization

Two different numerical grids are used to represent Ω0 and Ω1 , respectively. Grid points of Ω1 ∪ Γ are moved with the cell velocity field v1 . To advect the fluid domain the movement of the boundary grid points along Γ is harmonically extended onto Ω0 . The set of equations is split to subproblems on the different domains, which can be solved successively. An implicit Euler method is employed for the Navier-Stokes equations, an explicit Euler method is used to evolve B. In the following, time steps indices are denote by superscripts. The splitting scheme in time step m + 1 reads: 1. Solve the Navier-Stokes equations in the fluid (i = 0) with boundary condition v0m+1 = v1m ˜ (see step (6)) is subtracted in the advection operators in the on Γ. The grid velocity v same manner as the cell velocity vcell . 2. Update the left Cauchy-Green tensor B according to Eq. (5) by an explicit Euler time discretization. 3. Calculate λ1 , λ2 , κ, K and n from the position of boundary grid points along Γ, see Hu et m+1 δE m+1 δEbend m+1 stretch al. 40 for details. This allows computation of δEtension , and . δΓ δΓ δΓ 4. Solve the Navier-Stokes equations in the cell (i = 1) coupled to Eqs. (4),(5) with boundary condition Sm+1 1

· n = (1 −

ω)Sm 1

δEtension δEstretch δEbend · n + ω S0 · n − − − δΓ δΓ δΓ

m+1

on Γ.

The relaxation constant is set to ω = 0.25 to remove the stiffness of the coupling between flow and elasticity. Advection terms are neglected in this step since Ω1 is co-moving (see next step). 5. Move every grid point of Ω1 ∪ Γ with velocity v1 − vcell .

10


Page 11 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6. Solve the following Laplace problem to harmonically extend the interface movement into Ω0

∆˜ v=0

in Ω0 ,

˜=0 v

on δΩ0 \Γ, on Γ,

˜ = v1 − vcell v ˜. and move Ω1 including every grid point with velocity v

All equations are discretized in the corresponding axisymmetric formulation 45 such that they can effectively be solved in 2D. The meshes are created by gmsh, 46 the discretization with the Finite Element method is implemented in the software AMDiS. 47

4

Results

We simulate a typical RT-DC measurement of a single cell flowing through a square channel of width L = 20µm. To use the cortex stretch elasticity model given above, the flow is mapped onto an axisymmetric channel flow using the concept of the equivalent channel radius (see Sec 2). Assuming rotational symmetry of the cell deformation, axisymmetric differential operators can be used to effectively restrict the computations to two dimensions and largely reduce the computational complexity of the simulations. Hence, in the axisymmetric setting, a single spherical cell of radius r is initially placed in the middle of a circular channel of radius R = 0.547L. Since the channel moves with the cell barycenter, it suffices to consider a small portion of the channel length. We find that the cell deformation is not influenced by the length of the computational channel for a channel length of 40µm and larger. Accordingly, we specify Ω1 = {(x, y)|y ≥ 0, x2 + y 2 < r2 }, Ω0 = ([0, 40 µm] × [0, R]) \Ω1 . No slip boundary conditions are imposed on the channel wall. A pressure gradient pm is 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

imposed between channel inlet and outlet to drive the flow. Initially p0 = 2500 N/m2 , over time this pressure is adapted, pm+1 = pm /(0.9 + 0.1 · Q/(0.04 µl/s)), to yield a defined flow rate Q = 0.04 µl/s as in RT-DC experiments. The vertical velocity (in y-direction) is set to zero at the symmetry axis as well as at the channel inflow and outflow. The densities is set to ρ0 = ρ1 = 1000kg/m3 . The overall cell velocity is computed as vcell =

R

Ω1

v1 dV/

R

Ω1

1 dV. The time step is chosen

as dt = 0.2 µs (except for simulations for the hyperelastic shell model with E2D ≤ 0.05 N/m, where dt = 0.05 µs). The grid size is 0.2 µm around the cell surface and 1 µm in the bulk phases. As end time for the simulations we choose 500 µs except for the linear elastic bulk model with E = 15 kPa (1200 µs) and for the hyperelastic bulk model (2200 µs). The viscosity of the solvent is η0 = 0.015 Pa·s, as in Otto et al. 1 In RT-DC experiments, camera snapshots are taken at the end portion of the channel where the cells assume a stationary shape. Since the velocity gradient inside the cell vanishes in the stationary state, results are independent of the viscosity of the cytosol. Hence, for simplicity and numerical stability we choose η1 = η0 , throughout this work. Time-dependent simulations of channel entry and channel exit will depend on η1 . Such simulations may provide interesting additional information about the cell state and are left to future work.

4.1

Comparison with analytical model

While the numerical model can be used to combine bulk elasticity with cortex elasticity and surface tension, we will focus in the following on the two extreme cases that have also been investigated in Mietke et al.: 2 pure bulk elasticity and cortex elasticity with surface tension. In the latter case we fixed the surface tension to the elastic modulus of the shell, γ = 0.1E2D , to stick to a single free cell parameter. The bending rigidity is neglected for these first tests, cb = 0, and we will show in Sec. 4.3 that it has no influence on cell deformation in a reasonable parameter regime. Cell deformation is quantified by the imaged two-dimensional projection of the cell shape, as 12


Page 12 of 34

Page 13 of 34

it would be seen by a camera. Comparing the cell perimeter P with the perimeter of a circular cell of equal area A leads to the deformation measure d = 1 − 2

√

πA . P

Note that d = 0 for a circle

and d > 0 for any other shape. The time evolution of simulated cell deformation is depicted in Fig. 2(left). All four models, i.e. linear bulk elasticity, hyperelastic bulk material, linear shell elasticity and hyperelastic shell elasticity reach a stationary state of deformation within a few hundred microseconds. Note, that we do not simulate channel entry but start with an initially spherical cell within the channel. Hence, the given times only provide an indicator of the real time it takes for entering cells to reach a stationary state. In general, we observe that, in our simulations, the time until a stationary state is reached, correlates with cell stiffness as larger elastic moduli lead to earlier stationary states. Excellent volume conservation of simulated cells is shown in Fig. 2(right). 0.07

0.002

linear elastic shell E2D = 0.0313 N/m hyperelastic shell E2D = 0.0313 N/m linear elastic solid E = 3437 kPa hyperelastic solid E = 3437 kPa

0.05

linear elastic shell E2D = 0.0313 N/m hyperelastic shell E2D = 0.0313 N/m linear elastic solid E = 3437 kPa hyperelastic solid E = 3437 kPa

0.0015 relative volume change

0.06

deformation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.04 0.03 0.02 0.01

0.001 0.0005 0 −0.0005 −0.001 −0.0015

0

−0.002 0

100

200 300 Time [µs]

400

500

0

100

200 300 Time [µs]

400

500

Figure 2: Time evolution of deformation and volume of a single simulated cell of initial radius r = 7.66 µm for different models. Left: The cell deformation in the simulations assumes a stationary state after a few hundred microseconds. The four simulations displayed here are marked as black dots in Fig. 3 and Fig. 4. Right: The relative volume change is determined by (V − V0 )/V0 for volume V of the cell and initial volume V0 . The volume of the cell is very well conserved. Next, we compare stationary cell deformations of the analytical linear elastic model 2 with our numerical linear elastic model. As a typical output of RT-DC measurements, cell deformation, d, is plotted versus cell size (imaged cell area), where each cell of the measured population creates a single data points in a scatterplot. 1 So-called isoelasticity lines based on the analytical model 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

have been introduced to to identify regions of similar cell stiffness in such graphics. 2 Hence, it is very convenient to use these lines for a direct comparison of analytical and numerical model in the following. Fig. 3 shows the isoelasticity lines for both, elastic solid and elastic shell. The stationary deformations observed at the end time of Fig. 2 corresponds to single points in this plot, highlighted by filled black circles. The analytical model neglects nonlinear strain terms as well as the back-coupling of cell deformation on the flow field and therefore only makes reliable predictions for small deformations. Indeed, in the regime of small deformations d < 0.02 we find good agreement with our numerical model. For larger deformations and larger imaged cell areas (A ' 160µm2 ) deviations between the analytical predictions and our numerical model expectedly increase. This effect is even more pronounced for pure shell elasticity with surface tension where we find discrepancies between the analytical and the numerical model already for d ' 0.005.

4.2

Hyperelastic models

By including the back-coupling of cell deformation on the flow field, our numerical model naturally provides more accurate results for larger deformations. Due to the higher strains in such cases we employ neo-Hookean material laws for cell bulk and cortex. Such models are widely used to describe incompressible biological tissues, see Sec. 1. In Fig. 4 we compute isoelasticity lines for the neo-Hookean bulk and shell model and compare them to previous numerical results. For smaller deformations, d / 0.02, we find very good agreement between neo-Hookean and linear elasticity. For larger deformations neo-Hookean elasticity mostly leads to lower deformations than the linear elastic models. Generic cell shapes are depicted in Fig. 5. While linear elastic and neo-Hookean model yield very similar shapes, the analytical model shows larger deviations, in particular for more deformed cells. In contrast to predictions of the analytical model, 2 the numerical model does not show any lateral grooves for the elastic shell model. However, a dent occurs at the cell rear for larger deformations with the elastic shell model. The 3D side view in Fig. 5 (bottom left) depicts how the flatness of the trailing edge observed in experiments (cf. inset Fig. 9) could in principle arise 14


Page 14 of 34

Page 15 of 34

Elastic solid

0.06

E2D = 0.0836 N/m E2D = 0.0731 N/m E2D = 0.0627 N/m E2D = 0.0522 N/m E2D = 0.0418 N/m E2D = 0.0313 N/m E2D = 0.0209 N/m

0.1 0.08 deformation

0.08 deformation

Elastic shell

E = 11.45 kPa E = 9.16 kPa E = 6.87 kPa E = 5.73 kPa E = 4.58 kPa E = 3.44 kPa E = 2.29 kPa

0.1

0.04 0.02

0.06 0.04 0.02

0

0 60

80

100 120 140 160 180 200 220 240

60

80

area [µm2]

100 120 140 160 180 200 220 240 area [µm2]

0.02

0.02

0.015

0.015 deformation

deformation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.01

0.005

0.01

0.005

0

0 60

80

100 120 140 160 180 200 220 240

60

80

2

100 120 140 160 180 200 220 240 area [µm2]

area [µm ]

Figure 3: Isolines in area-deformation space for elastic solids (left) and elastic shells with a surface tension of γ = 0.1E2D (right). Lines mark the analytical results, 2 circles denote the results of the numerical linear elastic model. The two filled black circles in the top row indicate the solutions corresponding to Fig. 2. The bottom row shows a close-up of the respective top row diagram indicating good agreement between the models for small deformations.

15



Elastic solid

0.08

E = 11.45 kPa E = 9.16 kPa E = 6.87 kPa E = 5.73 kPa

Elastic shell

E = 4.58 kPa E = 3.44 kPa E = 2.29 kPa


0.1

deformation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 34

0.06 0.04 0.02


0.06 0.04 0.02

0

0 60

80

100 120 140 160 180 200 220 240

60

80

2

100 120 140 160 180 200 220 240 area [µm2]

area [µm ]

Figure 4: Isolines in area-deformation space for elastic solids (left) and elastic shells with a surface tension of γ = 0.1E2D (right). Crosses mark the results of the hyperelastic model, circles denote the results of the linear elastic model. The two filled black circles in the top row indicate the solutions corresponding to Fig. 2. from an effective convexification of indented cells due to the lateral imaging of cells in RT-DC. To make our results comparable with RT-DC experiments, we compute cell shape parameters, such as deformation and area, from cell contours that are convexified at the cell rear, throughout this work, see Fig. 5 (bottom right).

4.3

Influence of bending rigidity and finite Reynolds number

As pointed out above, the bending rigidity is believed to have very low impact on the cell deformation. To verify this hypothesis we conduct simulations with varying bending rigidity in the range cb = 0...2.5 × 105 kb T . We use the linear bulk elastic model with a fixed cell area of 115µm2 and fixed Young’s modulus, E = 2.29kPa. The resulting cell deformations are shown in Fig. 6 (left). The bending stiffness has no significant effect below cb ≈ 104 kb T . Since, typical values for lipid bilayer membranes are cb = 25...100kb T , this confirms that bending stiffness can be neglected. The Reynolds number Re = ρ0 vcell R/η0 is defined as the ratio of inertial forces to viscous forces and can be used to find a scaling invariance between two different cases of fluid flow. For the used parameter set, Re = 0.133, which is typically assumed small enough to render the 16


Page 17 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Elastic Solid

r = 4.38µm E = 9.16kPa

r = 7.11µm E = 9.16kPa

Elastic Shell

r = 5.69µm E = 2.29kPa

r = 5.47µm E2D = 0.0209N/m

3D Illustration of dent (in)visibility

r = 8.37µm E2D = 0.0209N/m

r = 7.11µm E2D = 0.0209N/m

r = 8.37µm E2D = 0.0209N/m

Convexified Linear Shell

r = 5.47µm E2D = 0.0209N/m

r = 7.11µm E2D = 0.0209N/m

r = 8.37µm E2D = 0.0209N/m

Figure 5: Top row: Cell shapes for elastic solid and elastic shell model for various parameters. Analytical solution 2 (blue circles), linear elastic model (red dots), neo-Hookean model (black line). Depicted cell sizes are matched for better visibility. All models agree for small deformations, while deviations from the analytical model become more prominent for larger deformations. Bottom row: The 3D illustration (left) shows a dent at the cell rear that is hidden in the lateral camera view. Accordingly, computed cell cross sections are convexified (right) for comparison with experiments, as shown here for the linear elastic shell model. inertial term in Eq. (1) negligible. 2 Our numerical model allows for an analysis of the effects that a finite Reynolds number has on cell deformations, which can be particularly useful if future RT-DC devices may use larger flow rates or smaller solvent viscosities. To test the effect of larger Reynolds number, we use the linear bulk elastic model with a fixed cell area of 115µm2 . The Reynolds number is changed by variation of solvent density in the range of ρ0 = 103 ...105 kg/m3 which yields an up to 100-fold increase in Re. As shown in Fig. 6 (right), we find a linear dependence of cell deformation on Reynolds number. If Re is increased by a factor of 10, the deformation increases approximately by only 1%. This holds for small deformations (E = 4.58kPa) as well as for large deformations (E = 2.29kPa). Similar results are obtained with the other elastic models. We conclude that inertial forces can very well be neglected in the original parameter setting and for a wide range of larger Reynolds numbers. In this case the Navier-Stokes equation (1) reduces to the Stokes equation, which is linear in the velocity. As a consequence, given that

17



0.055

0.064

0.054

0.062

E = 2.29 kPa 0.06

0.053

deformation

0.058 deformation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 34

0.052 0.051 0.05 0.049 0.048 0.047 100

1000 10000 cb [kbT]

100000

0.056 0.054 0.0126 0.0124 0.0122 0.012 0.0118 0.0116

E = 4.58 kPa

0

2

4

6

8

10

12

14

Re

Figure 6: The effect of bending stiffness and Reynolds number on deformation of a solid elastic cell. Left: The bending stiffness has no significant influence below cb ≈ 104 kb T . Right: A hundred-fold increase in Re changes the deformation by approximately only 10%. the channel is long enough to develop a Poisseuille flow profile at the inlet and outlet, the quasi-stationary state of the system is determined by only three non-dimensional parameters: ER3 /η0 Q, E2D R2 /η0 Q, γR2 /η0 Q. Hence, the universality of the isoelasticity lines found in the analytical model 2 carries over to the numerically computed isoelasticity lines, i.e., a change of experimental parameters only leads to rescaling of the involved elastic moduli without changing the lines in the plot. To be more precise, a change of (R, Q, η0 ) to (R′ , Q′ , η0′ ) can be accounted for by the replacing the parameters (E, E2D , γ) in the graphic by (ER/R′ , E2D , γ) × Q′ η0′ R2 /(Qη0 R′2 ) and scaling the cell area axis by a factor of (R′ /R)2 . This makes the isoelasticity lines a universal look-up graphic that can easily be adapted to a given set of experimental parameters. Plotting these lines together with data points from RT-DC experiments allows identification of regions of similar cell stiffness. 2

4.4

Comparison with experiments

We illustrate the potential to correctly predict elastic moduli by comparison of numerical results and RT-DC measurements. To this end, we use beads of polyacrylamide (PAAm) that can be considered as elastic solids. A sample of approximately 3000 beads was measured by RT-DC, see

18


Page 19 of 34

Fig. 7 (left). Note that an experimental flow rate Q = 0.08 µl/s and viscosity η0 = 0.0079 Pa·s has been used which leads to a different scaling of the iso-elasticity lines than before. The numerical hyperelastic solid model was employed with a bead size that corresponds to the mean of the experimental data (area ≈ 144 µm2 ). The numerical values were interpolated by a quadratic spline function providing a functional relationship between Young’s modulus and bead deformation. From this function, Young’s moduli were determined for beads of similar area (144 µm2 ± 2 µm2 ). The resulting distribution of predicted bead elasticity is depicted in Fig. 7(right). For comparison a subset of beads was measured by atomic-force microscopy (AFM). We find an excellent agreement between RT-DC and AFM measurements of bead elasticity.

0.06

20000 E= 2405 Pa E= 3607 Pa E= 4809 Pa E= 6011 Pa E= 7213 Pa E= 9618 Pa E= 12023 Pa Bead-Dat a

0.05

18000 16000 14000 E in Pa

0.04 deform at ion

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.03

0.02

12000 10000 8000 6000

0.01 4000 0.00 100

2000 120

140

160

180

200

RT-DC

AFM

Area in µm 2

Figure 7: Comparison of RT-DC and AFM measurements of PAAm beads. Left: The sample of beads was measured with RT-DC at a flow rate of 0.08 µls . Right: The boxplot shows excellent agreement of elastic moduli predicted by RT-DC and measured by AFM.

4.5

Determination of multiple cell parameters

Cell deformation as measured by the circularity of the cell, is very sensitive to the cell boundary contour extracted from camera images. Small deviations in the cell contour might occur by processing the noisy images and can increase the measured cell perimeter and lead to an overestimated deformation value. Integral quantities are typically more robust to such fluctuations. The moment-of-area tensor is such an integral quantity that can be used to quantify cell deformation 19



in the following. For a two-dimensional cell image, the moment-of-area tensor is given by 



2 (x − x0 )(y − y0 )   Ω (y − y0 ) Ω1 T= R 1 , R 2 (x − x0 )(y − y0 ) Ω1 (x − x0 ) Ω1

where x0 =

R

Ω1

R

R

x/|Ω1 | and y0 =

R

Ω1

(9)

y/|Ω1 | are the coordinates of the cell barycenter. The two

eigenvalues Λ1 , Λ2 of T are related to the length of the two semi-axes of the cell. The integration over the whole cell domain makes these values very robust with respect to uncertainties in the cell boundary contour. Note that, it is very easy to calculate Λ1 , Λ2 : If the cell image is horizontally R symmetric, the two eigenvalues are simply the diagonal elements of T: Λ1 = Ω1 (y − y0 )2 , Λ2 = R (x − x0 )2 . If the cell shape is given as a polygon of extracted contour points, the momentΩ1 of-area tensor can be easily evaluated as shown in the appendix Sec. 7.

Elasticity lines are shown for the neo-Hookean bulk and shell model for the ratio Λ1 /Λ2 in Fig. 8. For the solid elastic model we find no overlapping of lines for cells larger than 70µm2 , which confirms that this deformation measure is well-suited to conclude cell elasticity for a given deformation. For the elastic shell model we find significant overlapping of different isoelasticity lines, which makes it impossible to extract a definite E2D for a given Λ1 /Λ2 . 1.05

1.15

1 1.1 0.95 1.05

0.9 0.85 0.8 0.75 0.7 0.65 0.6

Λ1/Λ2

Λ1/Λ2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 34

E = 11.45 kPa E = 9.16 kPa E = 6.87 kPa E = 5.73 kPa E = 4.58 kPa E = 3.44 kPa E = 2.29 kPa

60

80

1 0.95 0.9 0.85 0.8

100 120 140 160 180 200 220 240


60

80

2

100 120 140 160 180 200 220 240 2

area [µm ]

area [µm ]

(a) Elastic solid

(b) Elastic shell

Figure 8: Isoelasticity lines for the new deformation measure Λ1 /Λ2 determined from the momentof-area tensor. However, the fact that Λ1 /Λ2 appears so different for the elastic shell and the elastic solid 20


Page 21 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


model can be very helpful to discriminate between the effects of these two models and to extract both parameters from a single cell. While we restricted our simulations so far to the two cases of either pure bulk elasticity or pure cortex elasticity with surface tension, the numerical model can be used to combine both effects. To illustrate this we conduct a parameter study where bulk and cortex elasticity are varied simultaneously. The cell size and surface tension are fixed to r0 = 7.27µm, γ = 0.1E2D . The results are plotted in a deformation space spanned by d and Λ1 /Λ2 in Fig. 9. An RT-DC image of an embryonic mouse stem cell of the same size is included and corresponds to a single location in this deformation space. Note, that the flatness of the trailing edge might indicate the presence of a hidden dent at the cell rear (cf. Fig. 5) which would also explain the rapid change in surface curvature at the back corners. By plotting lines for different bulk modulus E with colors indicating the shell modulus E2D , it is possible to estimate both elastic moduli of the experimental cells from nearby lines and colors. In our example the cell would have a bulk modulus of E ≈ 2.3kPa and a shell modulus of E2D ≈ 500pN/µm. In principle all mechanical parameters could be determined from shape data, which would require a more comprehensive sampling of the space of mechanical parameters and cell areas. This is not within the scope of this study, but currently addressed in a follow-up work.

5

Conclusion

In this paper, we presented numerical simulations of single cells in a circular channel that can be used to extract elastic cell parameters from RT-DC measurements. Therefore, we describe the cell as a visco-elastic bulk material enclosed by a thin cortical shell, subject to stretch elasticity, bending stiffness and surface tension. The cell bulk as well as the cell cortex are assumed as isotropic elastic materials, using either linear elasticity or neo-Hookean hyperelasticity. The discretization involves two separate Eulerian grids for the cell and the fluid domain that overlap at the cell surface. The full Navier-Stokes equations are solved within both domains, the corresponding surface stresses arising from shear and pressure forces enter as boundary conditions.

21



surface tension =0.1E h 0.1

E=763 Pa E=1146 Pa E=1909 Pa E=2291 Pa E=2673 Pa E=3055 Pa


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 34

0.07

0.008 0.007 0.006 0.005

0.06 0.05

0.004

0.04

0.003

0.03

0.002

0.02

0.001

0.01 0 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 E 2D [N/m] 1/ 2

Figure 9: A combined parameters study for r0 = 7.27µm, γ = 0.1E2D . The two deformation measures at the axis allow to discriminate between effects from bulk and cortex elasticity. Different lines correspond to different bulk moduli (in Pascal), while colors indicate the modulus of the cortical shell (in N/m). The inset shows a cell image acquired from an RT-DC measurement, as well as its position in the shape space (black dot) illustrating how mechanical properties can be determined for each cell.

22


Page 23 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The equations of motion are coupled to additional equations for the elastic stresses in the cell bulk and along the cell cortex by an explicit operator splitting. While the current simulations are restricted to axisymmetric scenarios, the model is easily extended to the fully three dimensional case. The only change would be the formulation of the membrane elasticity which can be included in 3D. 41 We applied the model to typical RT-DC measurements and thereby extend the first theoretical study on this topic that was based on an analytical model, which neglected the interaction between cell deformations and changes in the flow profile. 2 We used the concept of the equivalent channel radius 2 to obtain approximate cell deformations in a square channel. We could confirm that bending stiffness of the cortex/membrane complex has no effect on cell deformations for realistic values of bending rigidity. Therefore, we focused on the other mechanical elements of the cell and studied the influence of bulk elasticity and shell elasticity separately. We found that once the cell is in the channel, a stationary shape is reached within a few hundred microseconds, which provides an indicator of the total time between channel entry and stationarity. We compared the stationary cell shapes with the results of Mietke et al. 2 and found good agreement for small deformations where the analytical model is valid. For larger deformations (d ' 0.02) and larger cell areas (A ' 160µm2 ) the analytical model for bulk elasticity appears to leave its range of validity probably due to the missing back-coupling of cell deformation on the flow field. This effect was found even more pronounced for pure shell elasticity with surface tension where discrepancies between analytical and numerical model already occurred for d ' 0.005. In an experimental study with PAAm beads we illustrate the potential of the numerical method to predict elastic moduli from RT-DC measurements. We find an excellent agreement between the numerically predicted Young’s modulus and AFM measurements. By including the back-coupling of cell deformation on the flow field and keeping nonlinear strain terms, our numerical model naturally provides accurate results for larger deformations. Due to the higher strains in such cases we employed neo-Hookean material laws for cell bulk and cortex and provide isoelasticity lines that can be used for extraction of cell stiffness of significantly

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

deformed cells. We further found that inertial effects play no role in the original RT-DC parameter setting. This even holds for largely increased Reynolds number, which confirms the theoretical potential of RT-DC to work well with higher flow rates or less viscous carrier fluids. As a consequence the scaling invariance found in Mietke et al. 2 is carried over to the non-linear numerical model, i. e., a change of experimental parameters only leads to rescaling of the involved elastic moduli. Hence, the numerically computed isoelasticity lines are universal and can be scaled to account for any different flow rate, channel radius and viscosity. Therefore, using the concept of the equivalent channel radius our computed isoelasticity lines can be overlaid to data points of RTDC measurements to quickly identify regions of similar cell stiffness. Further, we introduced two new measures of cell deformation: the two eigenvalues of the moment-of-area tensor. Being integral quantities, these measures are very robust and are wellqualified to characterize deformation of even noisy cell images. We have shown that the quotient of these eigenvectors has the potential to distinguish between deformation effects stemming from cell cortex and cell bulk elasticity. While we restricted our simulations mostly to the two cases of either pure bulk elasticity or pure cortex elasticity with surface tension, the model can be used to combine all these effects. We demonstrated this here by using two different shape measures to simultaneously extract bulk elasticity and cortex elasticity from an experimental cell image. We are currently working on a combined study varying E, E2D and γ simultaneously. The produced cell shapes will be fitted to the three shape measures which should allow to extract the two elastic cell moduli and the cortical tension at once, and hence to disentangle the mutual contributions of these parameters to cell shapes in RT-DC. With this, the presented model provides an important stepping stone towards a quantitative, multivariate analysis of cell mechanical properties at highthroughput with many potential applications in medical diagnostics and biology. Acknowledgements SA acknowledges support from the German Science Foundation (grant AL 1705/3) and would like to thank the Isaac Newton Institute for Mathematical Sciences for its hospitality during the program Coupling Geometric PDEs with Physics for Cell Morphology,

24


Page 24 of 34

Page 25 of 34




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 34

Now calculate the time derivative of the shell stretching energy given in Eqs. (7)-(8),

dt Estretch =

Z

f ∇Γ · v + Γ

df ˙ df ˙ λ1 + λ2 dA, dλ1 dλ2

(11)

where the overdot denotes the material derivative and ∇Γ · the surface divergence operator. To p calculate λ˙1 , λ˙2 , we introduce the distance r to the symmetry axis x, r(x, y, z) = y 2 + z 2 and

the vector ¯r pointing away from the symmetry axis x, ¯r(x, y, z) = (0, y, z)/(y 2 + z 2 ). Note that the magnitude of ¯r scales like 1/r. From the definition of λ2 = r/r0 we conclude λ˙2 = r/r ˙ 0 = λ2 r/r ˙ = λ2 v · ¯r.

(12)

To obtain λ˙1 we consider the stretch of a surface area element A/A0 , which is related to the two principal stretches by A/A0 = λ1 λ2 . From mass conservation it is known that (A/A0 )˙ = A/A0 ∇Γ · v and hence (λ1 λ2 )˙ = λ1 λ2 ∇Γ · v

(13)

λ˙1 = λ1 (∇Γ · v − v · ¯r)

(14)

Substituting λ˙1 , λ˙2 in the energy time derivative yields

dt Estretch =

Z Γ

df f+ dλ1

∇Γ · v +

df df λ2 − λ1 ¯r · v dA, dλ2 dλ1

(15)

By splitting ∇Γ · v into ∇Γ · (Pv) + κn · v with projection operator P = I − n ⊗ n allows to use integration by parts on the ∇Γ · (Pv)-term to obtain

dt Estretch

df = v · (κn − ∇Γ ) f + dλ1 Γ Z

26

+

df df λ2 − λ1 dλ2 dλ1


· ¯r dA.

(16)

Page 27 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


With the definition of the energy densities in Eq. (8) we obtain δEstretch =(κn − ∇Γ ) δΓ  E2D    2(1−ν 2 )         =     E2D    3    

df f+ dλ1

+

df df λ2 − λ1 dλ2 dλ1

(17)

· ¯r

[(κn − ∇Γ )(3λ21 + λ22 + 4νλ1 λ2 + (1 + ν)(2 − 2λ2 − 4λ1 )) (thin shell)

− 2¯r(λ21 − λ22 + (1 + ν)(λ2 − λ1 ))]

−2 (κn − ∇Γ ) 3λ21 + 2λ22 + λ−2 r(λ21 − λ22 ) 1 λ2 − 6 − ¯

(neo-Hookean shell). (18)

In case of a simple thin shell we linearize the force to obtain δEstretch E2D [(κn − ∇Γ )(λ1 − 1 + νλ2 − ν) − ¯r(1 − ν)(λ1 − λ2 )] (linear shell). (19) = δΓ 1 − ν2

7

Calculation of moment-of-area tensor from a given cell contour

In experimental RT-DC data, the cell shape is usually described by a set of contour points that determine the cell area by a polygon. In the following, we show how the moment-of-area tensor can be calculated from these contour points without an explicit representation of the cell interior. Let (x0 , y0 ), ..., (xN −1 , yN −1 ) be the coordinates of the considered contour points, where N denotes the total number. We introduce the following values: • Area A of the cell polygon

A=

N −1 X

i=0 j=i+1 mod N

27

1 (xi yj − xj yi ) 2



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

• Integrals Ix =

R

x, Iy =

R

y over the polygon

1 Ix = 6 1 Iy = 6

N −1 X

(xi − yj − xj yi ) (xi + xj )

N −1 X

(xi − yj − xj yi ) (yi + yj )

i=0 j=i+1 mod N

i=0 j=i+1 mod N

• Coordinates xb , yb of the barycenter

xb =

Ix Iy , yb = A A

Now we define the point coordinates relative to the cell barycenter, (˜ xi , y˜i ) = (xi − xb , yi − yb ) for i = 0, ..., N − 1. The entries of the moment-of-area tensor can now be calculated by

T11

T22

1 = 12 1 = 12

T12 = T21 =

1 24

N −1 X

h

i i h x˜i y˜j − x˜j y˜i · y˜i2 + y˜i y˜j + y˜j2

N −1 X

h

i i h x˜i y˜j − x˜j y˜i · x˜2i + x˜i x˜j + x˜2j

N −1 X

h

i i h xi y˜i + x˜i y˜j + x˜j y˜i + 2˜ xj y˜j . x˜i y˜j − x˜j y˜i · 2˜

i=0 j=i+1 mod N

i=0 j=i+1 mod N

i=0 j=i+1 mod N

Note that for horizontally symmetric cell shapes, the eigenvalues of T are simply given by the diagonal elements: Λ1 = T11 , Λ2 = T22 .

References (1) Otto, O. et al. Real-time deformability cytometry: on-the-fly cell mechanical phenotyping. Nature methods 2015, 12, 199–202. (2) Mietke, A.; Otto, O.; Girardo, S.; Rosendahl, P.; Taubenberger, A.; Golfier, S.; Ulbricht, E.; 28


Page 28 of 34

Page 29 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Aland, S.; Guck, J.; Fischer-Friedrich, E. Extracting Cell Stiffness from Real-Time Deformability Cytometry: Theory and Experiment. Biophysical Journal 2015, 109, 2023–2036. (3) Ekpenyong, A. E.; Whyte, G.; Chalut, K.; Pagliara, S.; Lautenschläger, F.; Fiddler, C.; Paschke, S.; Keyser, U. F.; Chilvers, E. R.; Guck, J. Viscoelastic properties of differentiating blood cells are fate-and function-dependent. PLoS One 2012, 7, e45237. (4) Lautenschläger, F.; Paschke, S.; Schinkinger, S.; Bruel, A.; Beil, M.; Guck, J. The regulatory role of cell mechanics for migration of differentiating myeloid cells. Proceedings of the National Academy of Sciences 2009, 106, 15696–15701. (5) Henry, T.; Gossett, D. R.; Moon, Y. S.; Masaeli, M.; Sohsman, M.; Ying, Y.; Mislick, K.; Adams, R. P.; Rao, J.; Di Carlo, D. Quantitative diagnosis of malignant pleural effusions by single-cell mechanophenotyping. Science translational medicine 2013, 5, 212ra163– 212ra163. (6) Suresh, S. Biomechanics and biophysics of cancer cells. Acta Materialia 2007, 55, 3989– 4014. (7) Lincoln, B.; Erickson, H. M.; Schinkinger, S.; Wottawah, F.; Mitchell, D.; Ulvick, S.; Bilby, C.; Guck, J. Deformability-based flow cytometry. Cytometry Part A 2004, 59, 203– 209. (8) Guck, J.; Schinkinger, S.; Lincoln, B.; Wottawah, F.; Ebert, S.; Romeyke, M.; Lenz, D.; Erickson, H. M.; Ananthakrishnan, R.; Mitchell, D.; KÃďs, J.; Ulvick, S.; Bilbys, C. Optical deformability as an inherent cell marker for testing malignant transformation and metastatic competence. Biophysical journal 2005, 88, 3689–3698. (9) Remmerbach, T. W.; Wottawah, F.; Dietrich, J.; Lincoln, B.; Wittekind, C.; Guck, J. Oral cancer diagnosis by mechanical phenotyping. Cancer research 2009, 69, 1728–1732.

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(10) Hur, S. C.; Tse, H. T. K.; Di Carlo, D. Sheathless inertial cell ordering for extreme throughput flow cytometry. Lab on a Chip 2010, 10, 274–280. (11) Cooke, B. M.; Mohandas, N.; Coppel, R. L. The malaria-infected red blood cell: structural and functional changes. Advances in parasitology 2001, 50, 1–86. (12) Lee, G. Y.; Lim, C. T. Biomechanics approaches to studying human diseases. Trends in biotechnology 2007, 25, 111–118. (13) Swaminathan, V.; Mythreye, K.; O’Brien, E. T.; Berchuck, A.; Blobe, G. C.; Superfine, R. Mechanical stiffness grades metastatic potential in patient tumor cells and in cancer cell lines. Cancer research 2011, 71, 5075–5080. (14) Hou, H. W.; Bhagat, A. A. S.; Chong, A. G. L.; Mao, P.; Tan, K. S. W.; Han, J.; Lim, C. T. Deformability based cell marginationâĂŤa simple microfluidic design for malaria-infected erythrocyte separation. Lab on a Chip 2010, 10, 2605–2613. (15) Hou, H. W.; Li, Q.; Lee, G.; Kumar, A.; Ong, C.; Lim, C. T. Deformability study of breast cancer cells using microfluidics. Biomedical microdevices 2009, 11, 557–564. (16) Gossett, D. R.; Henry, T.; Lee, S. A.; Ying, Y.; Lindgren, A. G.; Yang, O. O.; Rao, J.; Clark, A. T.; Di Carlo, D. Hydrodynamic stretching of single cells for large population mechanical phenotyping. Proceedings of the National Academy of Sciences 2012, 109, 7630–7635. (17) Byun, S.; Son, S.; Amodei, D.; Cermak, N.; Shaw, J.; Kang, J. H.; Hecht, V. C.; Winslow, M. M.; Jacks, T.; Mallick, P.; Manalis, S. R. Characterizing deformability and surface friction of cancer cells. Proceedings of the National Academy of Sciences 2013, 110, 7580–7585. (18) Lange, J. R.; Steinwachs, J.; Kolb, T.; Lautscham, L. A.; Harder, I.; Whyte, G.; Fabry, B.

30


Page 30 of 34

Page 31 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Microconstriction arrays for high-throughput quantitative measurements of cell mechanical properties. Biophysical journal 2015, 109, 26–34. (19) Wyss, H. M. Cell Mechanics: Combining Speed with Precision. Biophys. J. 2015, 109, 1997–1998. (20) Dao, M.; Lim, C. T.; Suresh, S. Mechanics of the human red blood cell deformed by optical tweezers. Journal of the Mechanics and Physics of Solids 2003, 51, 2259–2280. (21) Eggleton, C. D.; Popel, A. S. Large deformation of red blood cell ghosts in a simple shear flow. Physics of Fluids (1994-present) 1998, 10, 1834–1845. (22) Zhou, E.; Lim, C.; Tan, K.; Quek, S. T. Finite element modeling of the micropipette aspiration of malaria-infected red blood cells. Microtechnologies for the New Millennium. 2005; pp 763–767. (23) Martins, P.; Natal Jorge, R.; Ferreira, A. A Comparative Study of Several Material Models for Prediction of Hyperelastic Properties: Application to Silicone-Rubber and Soft Tissues. Strain 2006, 42, 135–147. (24) Miller, K. Method of testing very soft biological tissues in compression. Journal of biomechanics 2005, 38, 153–158. (25) Lim H. W., G.; Wortis, M.; Mukhopadhyay, R. StomatocyteâĂŞdiscocyteâĂŞechinocyte sequence of the human red blood cell: Evidence for the bilayerâĂŞ couple hypothesis from membrane mechanics. PNAS 2002, 99, 16766–16769. (26) Zhang, J.; Johnson, P. C.; Popel, A. S. An immersed boundary lattice Boltzmann approach to simulate deformable liquid capsules and its application to microscopic blood flows. Phys. Bio. 2007, 4, 285. (27) Devendran, D.; Peskin, C. S. An immersed boundary energy-based method for incompressible viscoelasticity. J. Comp. Phys. 2012, 231, 4613–4642. 31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(28) Strychalski, W.; Copos, C. A.; Lewis, O. L.; Guy, R. D. A poroelastic immersed boundary method with applications to cell biology. J. Comp. Phys. 2015, 282, 77–97. (29) Fai, T. G.; Griffith, B. E.; Mori, Y.; Peskin, C. S. Immersed boundary method for variable viscosity and variable density problems using fast constant-coefficient linear solvers I: Numerical method and results. SIAM J. Sci. Comp. 2013, 35, B1132–B1161. (30) Peskin, C. S. The immersed boundary method. Acta numerica 2002, 11, 479–517. (31) Borazjani, I. Fluid–structure interaction, immersed boundary-finite element method simulations of bio-prosthetic heart valves. Comp. Meth. Appl. Mech. Engrg. 2013, 257, 103–116. (32) Dzwinel, W.; Boryczko, K.; Yuen, D. A. A discrete-particle model of blood dynamics in capillary vessels. J. Coll. Int. Sci. 2003, 258, 163–173. (33) Noguchi, H.; Gompper, G. Shape transitions of fluid vesicles and red blood cells in capillary flows. PNAS 2005, 102, 14159–14164. (34) Pozrikidis, C. Numerical simulation of the flow-induced deformation of red blood cells. Annals of Biomed. Engrg. 2003, 31, 1194–1205. (35) Sun, P.; Xu, J.; Zhang, L. Full Eulerian finite element method of a phase field model for fluidâĂŞstructure interaction problem. Computers & Fluids 2014, 90, 1–8. (36) Cottet, G.-H.; Maitre, E. A level set method for fluid-structure interactions with immersed surfaces. Math. Mod. Meth. Appl. Sci. 2006, 16, 415–438. (37) Lan, H.; Khismatullin, D. B. A numerical study of the lateral migration and deformation of drops and leukocytes in a rectangular microchannel. International Journal of Multiphase Flow 2012, 47, 73–84. (38) Turek, S.; Hron, J. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow ; Springer, 2006. 32


Page 32 of 34

Page 33 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(39) Salbreux, G.; Charras, G.; Paluch, E. Actin cortex mechanics and cellular morphogenesis. Trends in Cell Bio. 2012, 22, 536–545. (40) Hu, W.-F.; Kim, Y.; Lai, M.-C. An immersed boundary method for simulating the dynamics of three-dimensional axisymmetric vesicles in Navier–Stokes flows. J. Comput. Phys. 2014, 257, 670–686. (41) Barthes-Biesel, D.; Rallison, J. The time-dependent deformation of a capsule freely suspended in a linear shear flow. Journal of Fluid Mechanics 1981, 113, 251–267. (42) Zhu, Y.; Luo, X.; Ogden, R. W. Nonlinear axisymmetric deformations of an elastic tube under external pressure. European Journal of Mechanics-A/Solids 2010, 29, 216–229. (43) Höhn, S.; Honerkamp-Smith, A. R.; Haas, P. A.; Trong, P. K.; Goldstein, R. E. Dynamics of a Volvox embryo turning itself inside out. Physical review letters 2015, 114, 178101. (44) Huebscher, R. Friction equivalents for round, square and rectangular ducts. ASHVE Transactions (renamed ASHRAE Transactions) 1948, 54, 101–144. (45) Tomé, M.; Grossi, L.; Castelo, a.; Cuminato, J.; McKee, S.; Walters, K. Die-swell, splashing drop and a numerical technique for solving the Oldroyd B model for axisymmetric free surface flows. Journal of Non-Newtonian Fluid Mechanics 2007, 141, 148–166. (46) Geuzaine, C.; Remacle, J.-F. Gmsh: A 3-D finite element mesh generator with built-in preand post-processing facilities. International Journal for Numerical Methods in Engineering 2009, 79, 1309–1331. (47) Witkowski, T.; Ling, S.; Praetorius, S.; Voigt, A. Software concepts and numerical algorithms for a scalable adaptive parallel finite element method. Advances in Computational Mathematics 2015, 41, 1145–1177.

33




Page 34 of 34

Numerical Simulation of Real-Time Deformability Cytometry to Extract

Recommend Documents