2.4 Blood flow mechanics
2.4.3 Computational fluid dynamics
When systems become complex, maybe with an additional coupling to particle dynamics, an analytical solution of the Navier-Stokes equation (2.34) often cannot be achieved.
Therefore, fluid dynamics is often solved numerically and therefore a variety of methods in the field of computational fluid dynamics has evolved. In the following, the methods used in this thesis, the lattice-Boltzmann/immersed boundary and the boundary integral method, are shortly outlined.
Lattice-Boltzmann method
In contrast to other common methods such as finite volume or finite element method, the lattice-Boltzmann method does not discretize and solve the Navier-Stokes equation directly.
It is rather a mesoscopic approach which is based on the single particle distribution functionf(r,p, t) depending on position r and momentumpof a particle as well as time t [230–232]. The dynamics of the distribution function is determined by the Boltzmann equation. The Chapman-Enskog analysis proves that the Boltzmann equation leads to fluid behavior as governed by the Navier-Stokes equation [232, 233].
In the numerical realization, the lattice-Boltzmann method, a fluid is considered to be discretized on a regular Eulerian grid with nodes at positionxj and discretized velocities ci. We use the D3Q19 scheme with i= 0, . . . ,18 in three dimensions. The Boltzmann equation becomes the lattice-Boltzmann equation
fi(xj+ci∆t, t+ ∆t) =fi(xj, t) + Ωi(xj,Fj, t), (2.40) with Ωi the collision operator and ∆t a time increment. The collision operator Ωi governs changes in the distribution function due to collisions and further contains forces Fj, e.g., acting from a cell membrane onto the fluid. Those forces are transmitted to the fluid nodes by the immersed boundary method, as detailed below. We use a multiple relaxation time scheme [232] for the collision operator Ωi. Equation (2.40) can be integrated numerically including an update of the discrete distribution function due to the collision operator and propagation. The multiple relaxation time scheme further allows for the incorporation of thermal fluctuations of the fluid [234, 235]. Solid boundaries can easily be implemented by the bounce back scheme [232], whereas for elastic, deformable, and moving boundaries the lattice-Boltzmann method is accompanied by the immersed boundary method.
Immersed boundary method
The immersed boundary method provides a fluid-interface coupling for a deformable, moving boundary such as the cell membrane [236, 237]. The interface is represented
2 Modeling of biological membranes in the context of blood flow
by a Lagrangian grid of nodes immersed in the Eulerian fluid grid. In the present thesis, a discretization of the interface by nodes connected to triangles is used. From the corresponding constitutive equation the force per node f is calculated as, e.g., detailed in section 3.2 and in publication [pub1]. This force serves as input for fluid dynamics at the position of the membrane X0
Fj =Z f(X0, t)δ(X0−xj) d2X0, (2.41) and is mapped from the Lagrangian membrane mesh to the Eulerian fluid mesh by extrapolation from each membrane node to the surrounding fluid nodes. With the interface force as input, fluid dynamics is solved by a lattice-Boltzmann time step.
Afterwards, each interface node is advected with local fluid velocity, where the same interpolation scheme is used as for force spreading. All in all, the immersed boundary method provides a two-way, dynamic coupling of interface and fluid.
In the present thesis, the combined lattice-Boltzmann/immersed boundary method implemented in the software package ESPResSo [238, 239] has been used and extended.
Boundary integral method
In the limit of a small Reynolds number Re1, the boundary integral method can be used [110, 151, 240], which is directly based on the linearity of the Stokes equation (2.36) and as a consequence intrinsically neglects fluid inertia. Fluid dynamics is solved using the Green’s functions of the Stokes equation. Elastic interfaces such as the cell membrane are discretized and at each node on the discrete interface the traction jump in equation (2.24) calculated from the membrane forces is prescribed. Using the Green’s functions, the velocity in the simulation domain can be obtained for the current geometrical arrangement by the boundary integral equations. The interface nodes are then advected with the fluid velocity evaluated at the node position according to the kinematic boundary condition (2.39).
3 Numerical models for an active cell cortex
Following the thin shell theory from section 2.2, in section 3.1 the basics for the numerical representation of the cortex in a three-dimensional simulation is provided including the parabolic fitting procedure used to calculate geometrical properties and derivatives on the discrete cortex. In section 3.2, the computational method for an elastic, active cell cortex developed in publication [pub1] is introduced and in section 3.3 we summarize the method for a viscous cell cortex presented in publication [pub2].
σout
σin n
e1
e2
c
1 3 4 2
5
6 n
eξ eη
Fig. 3.1:Numerical membrane representation. Plasma membrane of the cell (blue) and cell cortex (red), as sketched on the left, are modeled together as thin shell (center) with the coordinate system (e1,e2,n) defined on the thin shell. Constitutive equations for the thin shell recover the elastic, viscous, or active physical properties of the membrane and the cortex. The thin shell is surrounded by an outer fluid and encloses an inner fluid described by the stress tensorσout andσin, respectively. In numerical simulations the thin shell is discretized by nodes connected to triangles (right). Considering its neighbors, a local coordinate system (eξ,eη,n) can be constructed on each node. Adapted from publication [pub1] with permission from APS and the left image from publication [pub5] with permission from the National Academy of Sciences..
3.1 Numerical representation of the cortex
In the following, we outline our numerical model for the cell cortex, which is based on a parabolic fitting procedure. As described in section 2.2, plasma membrane and cell cortex sketched on the left hand side of figure 3.1 are condensed into a deformable thin shell (center of figure 3.1). For its representation in numerical simulations we use a discretization of the thin shell by N points that are connected to triangles. This is a common approach for vesicle or cell simulations [69, 151, 240–242]. On the thin shell each node can be considered together with its neighborhood defined by the triangles adjacent to the central node. A node with its neighborhood is illustrated on the right hand side of figure 3.1.
As discussed above, thin shell theory is based on a coordinate system defined on the shell, e.g., in equation (2.3). The first step in the numerical procedure is the construction
3 Numerical models for an active cell cortex
of the local coordinate system at the position of each node rc
(eξ,eη,nc)rc. (3.1)
At first, the unit normal vectornc=n at the position of the central node is obtained by an average of the normal vectors nt on the adjacent triangles Tc
nc= X
t∈Tc
ntβt. (3.2)
The normal vector on each triangle nt is obtained by the cross product of the two vectors connecting the central node and its two neighbors. In equation (3.2) each normal vector is weighted by the angle βt of these two vectors. As the next step, for the first in-plane coordinate vector one neighbor node is chosen as reference, which is arbitrary but remains fixed during the simulation. The vector from the central node to the reference neighbor, e.g., x1, can be converted to the first in-plane unit vector by Gram-Schmidt orthogonalization and the second one can be obtained by a cross-product
eξ = x1−(x1·nc)nc
|x1−(x1·nc)nc|, eη = nc×eξ
|nc×eξ|. (3.3) Now, each vector from the central node can be expressed in the local coordinate system in equation (3.1) with in-plane coordinates ξ, η.
Following thin shell theory in section 2.2 we have so far obtained the coordinate system on the discrete shell and we can decompose an arbitrary vector along those coordinates.
This leaves us with the calculation of geometrical quantities, such as the metric tensor and the curvature tensor, but also the derivative of quantities defined on the shell. We achieve these using a parabolic fitting procedure such as used for calculation of bending forces [240]. The basic idea behind the parabolic fitting procedure is a Taylor expansion in the local coordinate system up to the second order around the central node. For an arbitrary function f the expansion in local in-plane coordinates (ξ, η) takes the form
f(ξ, η) =fc+Aξ+Bη+C
2ξ2+ D
2η2+Eξη, (3.4)
wherefc is the function value at the central node. Ato E represent the derivatives up to second order in the local coordinate system
A= ∇ξf|rc, B= ∇ηf|rc,
C = ∇ξ∇ξf|rc, D =∇η∇ηf|rc, E= ∇ξ∇ηf|rc, (3.5) which are evaluated at the position of the central node. Considering the function value fa at each of the N neighboring nodes and the theoretical expansion in equation (3.4) evaluated at the corresponding position (ξa, ηa) the squared residuum is defined as
χ2 = X
a∈N
[fa−f(ξa, ηa)]2 = X
a∈N
fa−fc−Aξa−Bηa− C
2ξa2− D
2η2a−Eξaηa
2
. (3.6) By minimization with respect to the parameter set (A, B, C, D, E)
∇(A,B,C,D,E)χ2 = 0, (3.7)