3.3 A numerical model for a viscous cell cortex
3.3.1 Velocity field on the discrete cortex
In the following, we consider a cortex in the long-time limit, where reorganization processes and motion of filaments under the action of motor proteins in the cortex take place [1].
On time scales larger than the viscoelastic relaxation time the active gel theory can be formulated in the viscous limit [20, 31, 37, 40, 162]. The total surface stress in equation
3.3 A numerical model for a viscous cell cortex
(2.22) thus consists of a viscous and an active stress
tα =tαvisc+tαact. (3.16)
The viscous surface stress consists of in-plane components tαβv
tαvisc=tαβv eβ, (3.17)
and we consider an in-plane, isotropic active stress in the following
tαact=ζgαβeβ =ζδαβeβ, (3.18) with its magnitude ζ =ζ(s1, s2) depending on the position on the membrane.
In contrast to the deformation field, which enters the method for an elastic cortex implicitly in the Skalak law (2.28) and Helfrich law (2.31), the key quantity here is the velocity field v of the cortex. The velocity field is a three-dimensional vector which can be expressed in local coordinates
v =vαeα+vnn with vα =v·eα and vn=v·n. (3.19) Analogous to the stress tensor in three-dimensional fluid dynamics in equation (2.33), in-plane gradients in the cortical velocity will contribute to the viscous stress. Therefore, we consider the in-plane velocity gradient
vαβ = 1
2[(∇αv)·eβ+ (∇βv)·eα]. (3.20) A key difference to three-dimensional fluid dynamics is the contribution of the normal velocity itself, rather than its gradient, to the velocity gradient in equation (3.20) by a coupling to the curvature tensor as can be seen from equation (2.19).
We take into account an in-plane shear viscosity ηs as well as an in-plane bulk viscosity ηb. Therefore, the constitutive law for the viscous surface stress becomes [9]
tαβv = 2ηs
vαβ − 1 2vγγgαβ
+ηbvγγgαβ. (3.21) Combining the viscous surface stress in equation (3.21) and the active surface stress in equation (3.18) according to equation (3.16), the force balance equations (2.25) and (2.26) for a force free cell cortex become
2ηs∇αvαβ + (ηb−ηs)∇αvγγgαβ =−∇αζαβ, (3.22)
−2ηscαβvαβ−(ηb−ηs)vγγcαβgαβ =−P +cαβζαβ, (3.23) where P is the pressure difference between inner and outer fluid. The force balance equations are accompanied by several constraints. We consider an isolated cell and as a consequence the total velocity and angular momentum are zero. In addition, the cytoplasm is incompressible. These three constraints are expressed mathematically by the integrals
Z
S
vdS=N−1X
ν=0 vνAν =0, Z
S
r×vdS=N−1X
ν=0 rν×vνAν =0,
Z
S
v·ndS =NX−1
ν=0
vν ·nAν =0,
(3.24)
3 Numerical models for an active cell cortex
where the first equality corresponds to the discretization of the cortex and Aν is the area per node.
This leaves us with the force balance equations valid on the whole membrane, i.e., for each node on the discrete membrane, and with the constraints. Together, they form a closed system of equations, which can be solved for the velocity field v and the pressure difference across the cortex P. The goal is to solve the viscous force balance equation for the velocity field v on the discrete cortex numerically. This is again based on a parabolic fitting procedure as detailed in section 3.1. Corresponding to equation (3.4) we expand the velocity vector up to second order around each node ν with velocity vν
v(ξ, η) =vν +∇ξv In contrast to before, the parabolic fitting procedure is not directly performed, because the actual velocity is now unknown. We rather invert the parabolic fitting procedure, i.e., the minimization of the theoretical expansion to the value at the neighboring nodes. We do so analytically to obtain the derivatives at node ν as a function of the velocity values of the node itself and its neighbors a(ν) = 1, . . . Nν
Av =Av
vν,{va(ν)}. (3.26)
The second essential step is to define a squared difference for the force balance equations (3.22) and (3.23) evaluated in the local coordinate system (cf. equation 3.1) at each node.
In the following, we abbreviate the left hand side as l.h.s.ρ and correspondingly the right hand side as r.h.s.ρ with the index ρ= 1,2, n enumerating the three force balance equations. Furthermore, the discretized constraints in equation (3.24) are considered with Lagrange multipliersλi and we define in total
χ2 =X
Here, the constraint of zero total velocity splitting up into its three spatial components is given for illustration. The χ2 is then expressed in terms of the velocity derivatives, which in turn depend on the actual velocity evaluated at all nodes. Together with the derivatives analytically depending on the velocity values according to equation (3.26), the χ2 can be minimized with respect to the velocity values
arg min
{vν}
χ2 = arg min
{vνx},{vνy},{vzν}
χ2. (3.28)
In stark contrast to the elastic method in section 3.2 where forces are locally and independently calculated, the coupled system of equations cannot be solved locally.
Rather the force balance equations, which are similar to a Poisson equation, have to be solved globally. The dimensions of the system of coupled equations excluding additional constraints is (3N + 1)×(3N+ 1) for the 3 components of the velocity vector evaluated at the N nodes plus the pressure.
3.3 A numerical model for a viscous cell cortex
This system of equations we write in matrix form as
∂χ2/∂v0x
Numerically solving this system of equations (3.29), e.g., by LU decomposition, gives the velocity field v on the cortex as set of the velocity values at the nodes {vν}.
3.3.2 Validation
In order to validate the numerical model of the cell cortex in the viscous limit we choose two different approaches: in the static case the velocity field obtained on a spherical cortex is compared to analytical solutions, which we obtain for an active stress distribution in terms of spherical harmonics. This further allows for a detailed comparison of the different quantities in simulations with the theoretical results. Next, we consider a dynamically deforming cell cortex with an active stress increasing around the equator, where we can compare the three-dimensional dynamic deformation and evolving flow field on the cortex to axisymmetric simulations [249, 250] performed by our collaboration partners in publication [pub2].
First, we apply an active stress expressed in spherical harmonics ζ(θ, φ) =ζ0+X
n,m
ζnmYnm(θ, φ), (3.30)
with a constant offset active stress ζ0, Ynm being the spherical harmonics and ζnm the expansion coefficients. Expanding the velocity field on the cortex in terms of vector spherical harmonics, expressing the viscous surface stress tensor as function of tensorial spherical harmonics and applying identities of those, we are able to derive analytical conditions to the expansion coefficients from the viscous force balance equations. The analytical conditions then give an analytical expression for the velocity field for given active stress distribution.
The analytical procedure itself can be used for two different test setups: the full force balance equations which lead to a pure normal velocity on the cortex is tested on the one hand. On the other hand, the normal force balance can be replaced by a condition fixing the normal velocity to zero. The latter, somewhat artificial test setup allows us to obtain an analytical solution for a finite tangential velocity along the cortex. In both scenarios we can apply the algorithm described above and obtain numerical results in excellent agreement with the analytical solution. This is done for an axisymmetric as well as a non-axisymmetric active stress in terms of Y20 and Y21, respectively. The three-dimensional velocity profile on the spherical cortex for the full system in case of the axisymmetric active stress in terms of Y20 with ηs = 1 and ηb = 1 in simulation units is
3 Numerical models for an active cell cortex
Fig. 3.3:Validation for a viscous active cortex. a) Three-dimensional velocity field numerically obtained on a static spherical cortex subjected to an axisymmetric active stress distribution in terms of the spherical harmonicsY20with the arrows giving the velocity direction and the color coding for the velocity magnitude. b) The numerically obtained velocity with pure normal component depending on the polar angleθis in very good agreement with the corresponding analytical solution. c) Initial and d) final velocity field on an evolving cortex with viscosities ηs= 1 andηb= 1 subjected to an active stress distribution according to equation (3.31) with exponent p = 4, widthσ = 10, offset ζ0 = 1, and magnitude ˆζ = 1. e) The dynamics in terms of the pole to pole and equatorial furrow radius (dots) agrees very well with results from axisymmetric simulations (lines). Time is considered relative to the time scaleta= ηζs0. Reprinted from publication [pub2].
3.3 A numerical model for a viscous cell cortex
shown in figure 3.3 a). The comparison of the numerically obtained normal velocity with the analytical solution over the polar angle is shown in figure b), where we obtain a very good agreement. These scenarios further show that the presented algorithm converges properly with the error decreasing systematically with increasing resolution.
In a second step, we consider an isotropic, axisymmetric active stress which increases around the equator according to
ζ(θ, φ) = ζ0+ ˆζexp−σ
θ− π 2
p
, (3.31)
with the exponent p being an even number, the magnitude ˆζ and the widthσ. For given active stress, we solve for the velocity field on the discrete cortex and use the normal velocity for an update of the cortex shape by applying the Euler algorithm. For an active stress distribution with exponent p= 4, offset active stress ζ0 = 1, and amplitude ˆζ = 1 with the viscosities ηs = 1 andηb = 1 in simulation units, we show the velocity profile on the initial and the finally deformed cortex in figure 3.3 c) and d), respectively. The prescribed active stress leads to a contraction around the equator and, corresponding to the incompressibility of the cortex interior, to an extension at the poles. The dynamics of both the equatorial contraction and the pole extension is in very good agreement with results from the axisymmetric simulations of our collaborators as shown in figure 3.3 e).
In addition, we compare the velocity field on the deformed, three-dimensional cortex to the one obtained by axisymmetric simulations at different times. Both tangential and normal velocity show very good agreement. In total, this successfully validates our three-dimensional model of a cell cortex in the viscous limit.
As a first test application of a deforming cell cortex in a non-axisymmetric situation, we apply our developed method to consider an initially spherical cortex with an active stress according to equation (3.31) and subjected to an initial shear deformation. In ref.
[pub2], we analyze the cortex evolution and the velocity field on the deforming cortex over time.
4 Rayleigh-Plateau instability for anisotropic tension and blood platelet formation
The Rayleigh-Plateau instability is a well known phenomenon, which describes the instability leading to break-up of a liquid jet into droplets [54]. It is known to be driven by the surface tensionγ of the jet interface [54, 251]. The surface tension leads to a pressure difference between the inner and outer fluid of the jet according to the Young-Laplace equation [54,150]
p=γ 1 Rφ
+ 1 Rz
!
, (4.1)
whereRφandRz are the radius of curvature in azimuthal and axial direction, respectively.
Initially, the jet is cylindrical with an unperturbed radius R0. A periodic undulation of the jet interface with wavelength λ along its axis leads to a spatial variation of the pressure p. Considering the first term of the Young-Laplace equation (4.1), at the site of a neck, the pressure contribution Rγ
φ increases due to a decrease of the local radius.
Correspondingly, Rγφ decreases at the site of a bulge. This mechanism amplifies the undulation even further. In contrast, the second term Rγz leads to an opposing effect: the curvature in axial direction leads to a positive pressure perturbation at a bulge and a negative one at a neck. This stabilizes the interface. The balance of both effects leads to a dominant perturbation mode, at which the interface becomes unstable and which grows fastest. For an ideal fluid jet with passive ambient fluid, the fastest growing wave number km= λ2πm fulfills the condition [51, 54]
kmR0 = 0.697. (4.2)
Based on the work of Plateau [50], a dispersion relation was first derived by Rayleigh [51], which yields the criterion in equation (4.2). For a viscous liquid jet in air again Rayleigh [252] and Chandrasekhar [253] derived a dispersion relation, which has been generalized to include an external fluid of arbitrary density and viscosity by Tomotika [52]. Experiments [55] have first shown that a Rayleigh-Plateau instability also occurs for tubular vesicles under externally induced tension [55–58,164, 254]. Considering cell or tissue tubes, the active stress in the cortex has been proposed to trigger a Rayleigh-Plateau like instability [32, 33], which we confirmed by three-dimensional simulations [pub1]. Indeed, an active stress in the cell cortex bares a striking similarity to surface tension as it enters the in-plane component of the surface stress in equation (2.22) in the same way (cf. equation (2.32)) as surface tension does in case of a liquid jet
tactαβ = Ta 0 0 Ta
!
, tjetαβ γ 0 0 γ
!
, (4.3)
4 Rayleigh-Plateau instability for anisotropic tension and blood platelet formation
for isotropic, contractile active stress Ta > 0. Therefore, despite their fundamental difference both isotropic active stress and surface tension lead to the same physical behavior. This analogy we use in the next section, where we refer to the term interfacial tension in general. Interestingly, there is a remarkable difference between cortical stress in cells or tissues and the classical scenario of a liquid jet with isotropic surface tension:
cells often have an anisotropic, cortical active stress [20, 31, 40,59–65].
In section 4.1, we summarize the effects which a stress anisotropy has on the Rayleigh-Plateau instability as investigated in detail in publications [pub3] and [pub4]. We combine analytical linear stability analyses with computer simulations based on our developed method [pub1]. The break-up of a liquid jet has a striking similarity to cell shapes occurring during blood platelet formation. Using our computational model [pub1], in section 4.2, we show that a biological Rayleigh-Plateau instability can explain the flow accelerated formation of blood platelets [pub5].