Image-based model of the spectrin cytoskeleton for red blood cell simulation

We simulate deformable red blood cells in the microcirculation using the immersed boundary method with a cytoskeletal model that incorporates structural details revealed by tomographic images. The elasticity of red blood cells is known to be supplied by both their lipid bilayer membranes, which resist bending and local changes in area, and their cytoskeletons, which resist in-plane shear. The cytoskeleton consists of spectrin tetramers that are tethered to the lipid bilayer by ankyrin and by actin-based junctional complexes. We model the cytoskeleton as a random geometric graph, with nodes corresponding to junctional complexes and with edges corresponding to spectrin tetramers such that the edge lengths are given by the end-to-end distances between nodes. The statistical properties of this graph are based on distributions gathered from three-dimensional tomographic images of the cytoskeleton by a segmentation algorithm. We show that the elastic response of our model cytoskeleton, in which the spectrin polymers are treated as entropic springs, is in good agreement with the experimentally measured shear modulus. By simulating red blood cells in flow with the immersed boundary method, we compare this discrete cytoskeletal model to an existing continuum model and predict the extent to which dynamic spectrin network connectivity can protect against failure in the case of a red cell subjected to an applied strain. The methods presented here could form the basis of disease- and patient-specific computational studies of hereditary diseases affecting the red cell cytoskeleton.

with the corresponding free energy given by where C is a constant independent of r. The force is given by F = −∇A, which gives i.e. a Hookean spring with zero rest length and spring constant k = 3k B T /(2Ll p ). Since the source of the elasticity is purely combinatorial, i.e. there are many more polymer conformations in which the end-to-end length is small rather than large, this is known as an entropic spring model. Although (3) displays reasonable behavior at small extensions [1], it breaks down at extensions approaching the contour length. Whereas the end-to-end distance of a polymer consisting of finitely-many monomers is at most N l k , the linear approximation above allows in principle for infinite extensibility. The phenomenological worm-like chain model corrects this by introducing a divergent force at finite extension [5]. To represent each spectrin polymer we therefore use the worm-like chain model, which has the force-extension relationship F wlc (r) = F wlc ( r ) r r , F wlc ( r ) = k B T l p The one case in which we use F linear is for the test of red cells with dynamic networks subject to prescribed strain in the section titled "Simulating response to flow and applied strain". When using the worm-like chain model, the number of irreversibly broken edges is difficult to detect above the noise of network dynamics, since the divergent force prevents edges from approaching the contour length (see the section below titled "Dynamic connections"). Note that F wlc agrees with F linear in the limit of small extensions, and arbitrarily large extensions do not occur when there is irreversible breakage beyond a threshold.

Density through slab
Given a slab with thickness η, we wish to compute the vertically-averaged node distance PDF (D), defined in the section "Generating the cytoskeleton" of the Main Text by b a (D)dD = Pr(D ij ∈ (a, b)), where D ij is the distance between any two nodes i and j.
As before, the nodes are assumed to be uniformly distributed through the slab. Throughout this appendix, we ignore the effects of edges in the horizontal directions. As η → ∞ we should have (D) → 4πD 2 /V , where V is the volume of the slab. However, we are interested in the case of a thin slab with η ∞, so the effects of finite η must be accounted for. We will consider several cases. First, take D < η/2, so that there is at most one intersection between a sphere of radius D centered inside the slab and the slab boundary. If the sphere is centered at a height h > D above the lower slab boundary, it is fully contained within the slab. If h < D, the sphere intersects the lower boundary at some angle θ (see Fig A). The excluded solid angle Ω is given by which implies that the excluded surface area is 2πD 2 (1 − cos θ). The vertically-averaged density is therefore given by where we have used the symmetry of the top and bottom halves of the slab. Since cos θ = h/D, Next, take η/2 <= D < η so that the sphere may intersect with at least one and possibly both slab boundaries. Assume without loss of generality that the sphere intersects the lower boundary. Let θ be as above and, in the case that the sphere intersects with the upper boundary, define φ to be the associated angle of intersection (see Fig B). Now, we must account for potentially two excluded regions of the sphere, which yields Note that the expression above is identical to the formula obtained for D < η/2. Finally, suppose that D > η. In this case, there is always an intersection with both upper and lower boundaries for all h and In summary, (D) is a piecewise function satisfying

S1 Text Dynamic connections
Here, we extend our model to allow for dynamic connections that change over time.
The dynamics of network connections are governed by rate constants k on (D) and k off (D) of formation and breakage, respectively. We assume that the rate constant of breaking is independent of distance so that k off (D) ≡ k off . We next show how to choose the function k on (D) that results in a steady-state distribution of edge lengths. Recall from the section "Generating the cytoskeleton" of the Main Text that p(D) is the probability that a pair of nodes separated by a distance D is connected by an edge. Now letting p = p(D, t) also depend on time, the dynamics of p for a cell at rest follow If we specify a constant rate of breaking k off (D) ≡ k off , we may use the steady-state condition and the form of p(D), derived in the section "Generating the cytoskeleton" of the Main Text, to set k on (D). Note that (12) is only used to derive the rate constant k on (D) of edge formation and we assume p(D) does not depend on time in all our simulations.
The assumption that k off is independent of D is arbitrary-the remodeling of spectrin is still not well-understood and there is insufficient experimental data to discriminate between competing models of network dynamics. However, for the purpose of this work the rate of remodeling matters more than its specific form and we believe the same qualitative effects would result for other models of network dynamics as well.
In our model the network dynamics are handled in a stochastic manner. To check whether new bonds are formed, for each eligible pair of points (i.e. not already connected and with separation D < D max ) we generate a random number uniformly between 0 and 1 at each time step. We draw a new edge between points if the corresponding random number is less that k on (D)∆t. To determine whether two nodes become disconnected, we generate random numbers for each edge and remove that edge if the random number is less that k off ∆t. For computational efficiency, we recycle the binning of nodes into boxes used for generating the model cytoskeleton in the section "Generating the model cytoskeleton on a triangulated surface" of the Main Text. That is, at each time step, we update the boxes containing the nodes using methods described in [6]. This leads to a considerable computational savings, since when making new connections only pairs of nodes within the same and adjacent boxes must be checked.
To characterize the viscoelastic behavior of this dynamic network we simulate two tests of viscoelasticity on a remodeling network of worm-like chains. This is done on a 2D patch of cytoskeleton. The network is expected to behave like a dashpot representing the viscosity of the ambient fluid in parallel with a Maxwell element, i.e. a spring in series with another dashpot. The spring and dashpot components of the Maxwell element represent the elastic polymers and the viscous effect of remodeling, respectively. This is referred to as a Fluid 1 element in [7], in which the authors systematically study the response of 3D cytoskeletal networks made up of different viscoelastic elements.
Following the descriptions of viscometer tests in [7,8], we first simulate a stress relaxation test in a viscometer by shearing the network and observing how the force required to maintain the shear decays over time (Fig C). To do this, a 2D random network made up of N = 1, 200 nodes along with its corresponding triangulation are generated as in the section "Elastic response of random graphs" of the Main Text and then subjected to a prescribed shear strain. Nodes located in layers on the top and bottom boundaries are fixed in place, while the interior nodes move according to overdamped dynamics determined by worm-like chain and area-preserving forces. To separate the timescale of the fluid from the timescale of remodeling, the network is fixed with a static connectivity and relaxed to equilibrium before any bonds are allowed to break and reform. From the ratio of the equilibrium force and displacement, the material has an effective elastic modulus of k ≈ 0.02 dyn/cm. Once network dynamics are turned on, the force decays with characteristic timescale of ≈ 0.04 s (Fig C(a)). Since this timescale is given by η 1 /k for a Fluid 1 element, where η 1 is the dashpot viscosity in the Maxwell element [7], i.e. the effective viscosity of network we find η 1 ≈ 8 · 10 −4 g/s for the parameters used. Next, we simulate a creep test in which, rather than a prescribed shear strain, we subject the network to a constant external force. The network is prepared as above. Nodes in a thin strip on the bottom boundary are held in place, but instead of prescribing the positions of the nodes on the top boundary, they are subjected to a constant force. In the case of a purely elastic material, the top nodes would not move past their equilibrium positions. This is indeed observed when the network is not allowed to remodel (Fig C(b)). However, with remodeling, there is a plastic deformation, or creep, and the nodes drift gradually in the direction of the applied force. The rate of drift is equal to (η 1 + η 2 ) −1 , the inverse of the sum of dashpot viscosities in the Fluid 1 element [7]. After the initial transient, we find η 1 + η 2 = 1.2 − 2.0 · 10 −3 g/s.
Putting these calculations together gives the estimate η 2 = 0.4 − 1.2 · 10 −3 g/s for the viscosity on the whole structure due to the fluid. This result is self-consistent, since the net friction of 1, 200 · 3.3 · 10 −7 = 4.0 · 10 −4 g/s, calculated by multiplying the friction coefficient ζ = 3.3 · 10 −7 times the total number of nodes, lies within this range.
To justify using a linear polymer force model instead of a worm-like chain for the optical tweezer simulations in the section "Simulating response to flow and applied strain" of the Main Text, here we compare differences in irreversible breakage events on a 2D patch of cytoskeleton using both the worm-like chain and linear spring models. The two ends of the domain are pulled apart at a fixed rate, which stretches out the material in a manner reminiscent of optical tweezers experiments. We consider a patch of cytoskeleton on a 2D inextensible sheet being stretched out in the presence of network dynamics for both the linear spring and worm-like chain models. Here we use a threshold of 0.9L for irreversible breakage. Upon simulation, we find many more broken edges when using the linear spring model, which is expected given the stiffening of the worm-like chain model at large extensions. Although the edge-length distribution extracted by the image processing has a small tail beyond the irreversible breakage PLOS 5/9 S1 Text

Immersed boundary method
The immersed boundary formulation is a mathematical description of fluid-structure interaction that combines the Navier-Stokes equations with the dynamics of an immersed elastic material (see for instance [9,10] and references therein). Fluid with density ρ and viscosity ν is assumed to fill a domain, and the fluid is described in cartesian coordinates x by a velocity field u(x, t) and pressure field p(x). The domain also contains an elastic material, which is convenient to describe in Lagrangian coordinates q that denote material points. There is a function X(q, t) that returns the Cartesian position of material point q at time t, and the elastic energy of the material can be computed as a functional of the deformed configuration so that E = E[X(·, t)].
This gives rise to a force density F(q, t) = −δE/δX, where δ/δX denotes the variational derivative with respect to the configuration of the elastic material.

S1 Text
The equations that characterize the immersed boundary formulation are as follows: f Equations (13)- (14) are the Navier-Stokes equations and (17) gives the force density on the membrane in Lagrangian coordinates. Equations (15)- (16) are the interactions equations that couple the fluid and elastic material. The Lagrangian force density is convolved against the Dirac delta function δ(x), giving rise to the force density f (x, t) in cartesian coordinates that appears in (13). In our case, the Lagrangian force involves one contribution from the membrane and another from the cytoskeleton. At a discrete point on the cytoskeleton, the force is computed by where the sum is over all points X j attached by an edge in the network to X k . At a discrete point on the membrane, on the other hand, the force is defined in terms of the energy W : where W = W bulk + W bend and ∇ X k denotes the gradient of the energy with respect to the position X k . We remark that F k is the force, not the force density, applied at point X k to the fluid. Finally, to make the formulation complete, we specify the membrane energy via where Z is the reference configuration, H is the total curvature, i.e. the sum of the two principal curvatures, and k bend = 2.0 × 10 −19 J [11] is the bending rigidity. The bulk modulus κ and Jacobian J are as in the section titled "Elastic response of random graphs". The total force is obtained by summing over all points on the membrane and cytoskeleton. See [10,12] for details regarding how these forces are computed in practice. The elastic material moves at the local fluid velocity according to (16). The immersed boundary equations describe the nature of the coupling between the membrane and cytoskeleton, which are treated as independent elastic structures in our simulations that are both immersed in the same ambient fluid. The forces they exert are both spread onto cartesian coordinates so that f (x, t) is the resultant force, which determines the fluid velocity in which both structures move. Because the membrane and cytoskeleton initially lie on the same surface (see Fig E), they stay together since they are advected in the same velocity field. Diffusion of cytoskeletal proteins and slip relative to the membrane can be included within the immersed boundary framework [13,14], and adding these effects would make the model even more realistic and allow for simulating an even larger range of observed phenomena, such as the nontrivial interplay between diffusion of cytoskeletal proteins and disruptions in both S1 Text vertical connections to the membrane and horizontal cytoskeletal connections observed in computational [15] and experimental [16] studies.
In the present context of the continuous problem, writing the interpolation equation as above with convolution against a Dirac delta function is superfluous; the velocity field u needs only to be evaluated at a material point X(q, t) to find the velocity there. However, this equation becomes useful once the equations are discretized, since replacing the Dirac delta function by a regularization yields a numerical method for performing the interpolation. The choice of regularized delta function is the key ingredient of the immersed boundary method; it should interpolate smooth functions accurately and have finite support for computational efficiency. In the simulations presented here we have used a delta function with a support of 4 gridpoints in each spatial dimension. For further details, see [9,10,17].
In our simulations, we have taken water to be the ambient fluid, so that ρ = 1 g/cm 3 and ν = 1.2 cP, and the domain to be a rectangular box of size 10 µm ×10 µm ×20 µm with periodic boundary conditions in each direction. A discretization of the box with either 128 × 128 × 256 or 256 × 256 × 512 is used, which results in an equal grid spacing ∆x in each direction. We have found it necessary to take a timestep ∆t with ∆t = (0.01 − 0.1)∆x to ensure stability. A predictor-corrector timestepping scheme is used [10], and the linear system resulting from discretizing the Navier-Stokes equations are solved with fast Fourier transforms. The surface of the red cell is represented by a triangulated surface composed of approximately 20, 000 triangles. The method has been implemented with shared memory parallelism on the GPU using Python and OpenCL. Typical simulations of cells with static networks take 1-7 days of continuous computing time depending on the resolution of the fluid grid, whereas typical simulations of cells with dynamic networks take several days to weeks.