Continuum modeling of non-conservative fluid membrane for simulating long-term cell dynamics

Abstract

Living cells actively deform and move by their force generations in three-dimensional (3D) space. These 3D cell dynamics occur over a long-term time scale, ranging from tens of minutes to days. On such a time scale, turnover of cell membrane constituents due to endocytosis and exocytosis cannot be ignored, i.e., the surface membrane dynamically deforms without mass conservation. Although membrane turnover is essential for large deformation of cells, there is no computational framework yet to simulate long-term cell dynamics with a non-conservative fluidic membrane. In this paper, we proposed a computational framework for simulating the long-term dynamics of a cell membrane in 3D space. For this purpose, in the proposed framework, the cell surface membrane is treated as a viscous fluid membrane without mass conservation. Cell shape is discretized by a triangular mesh, and its dynamics are expressed by effective energy and dissipation function. The mesh structure, distorted by membrane motion, is dynamically optimized by introducing a modified dynamic remeshing method. To validate the proposed framework, numerical simulations were performed, showing that the membrane flow is reproduced in a physically consistent manner and that the artificial effects of the remeshing method were negligible. To further demonstrate the applicability of the proposed framework, numerical simulations of cell migration induced by a mechanism similar to the Marangoni effect, i.e., the polarized surface tension actively generated by the cell, were performed. The observed cell behaviors agreed with existing analytical solutions, indicating that the proposed computational framework can quantitatively reproduce long-term active cell dynamics with membrane turnover. Based on the simple description of cell membrane dynamics, this framework provides a useful basis for analyzing various cell shaping and movement.

Graphical abstract

Appendix

1.1 A1: Membrane flow under area constraint

The motion equation of the vertices in Eq. (7) is derived from that of the continuum in Eq. (1); however, whether membrane dynamics can be reproduced in a physically consistent manner remains to be determined. Therefore, whether planar membrane flow can be reproduced when the total membrane surface was conserved was determined. Numerical simulations of planar membrane flow were performed, and the obtained velocity profiles were evaluated and compared with theoretical solutions.

In $xyz$-coordinates, a rectangular system was considered whose lengths along the individual axes are represented by $L_{x}$, $L_{y}$, and $L_{z}$ (Fig. 

Fig. 8

Hagen–Poiseuille flow of membrane. a. Snapshot of membrane flow resulting from numerical simulations. Triangles outside boundaries are not shown. Colors represent membrane velocity along the $x$-axis, $v_{x}$. b. Spatial profile of membrane velocity along the $x$-axis, $v_{x}$, under the condition with $F = 0.01\,\mu l^{2} /\tau$. Black dots indicate individual vertex velocities at $t = 100\tau$. The red curve represents the theoretical solution of Eq. (26). c–e. Average maximum velocity, $v_{{{\text{max}}}}$, as a function of pressure gradient, discretization length $l_{{\text{d}}}$, and threshold ratio $\lambda$, respectively. The red line in c-e represents the theoretical solutions of Eq. (27). Error bars in c-e indicate standard deviations

8a). A periodic boundary condition was imposed on $x = 0, L_{x}$ and $y = 0, L_{y}$. The membrane spreads in the $xy$-plane across the boundaries. By representing the $x$-component of the ith vertex location by $x_{i}$, the energetic force is introduced by the effective energy $U$ in Eq. (7), which is expressed as follows:

$$ U \equiv \left\{ {\begin{array}{*{20}c} { - \mathop \sum \limits_{i}^{{{\text{vertex}}}} \frac{F}{{\rho^{*} }}x_{i} } & {y_{i} < \frac{{L_{y} }}{2}} \\ {\mathop \sum \limits_{i}^{{{\text{vertex}}}} \frac{F}{{\rho^{*} }}x_{i} } & {y_{i} \ge \frac{{L_{y} }}{2}} \\ \end{array} } \right., $$

(24)

where constant $F$ is the stress exerted on the membrane that induces the pressure gradient as

$$ \frac{dp}{{dx_{i} }} = - \rho_{*} \nabla_{i} U. $$

(25)

Theoretically, this pressure gradient induces a Hagen–Poiseuille flow, whose velocity profile is given as

$$ v_{x} \left( y \right) = \left\{ {\begin{array}{*{20}c} { - 16v_{{{\text{max}}}} \frac{{y_{i} }}{{L_{y} }}\left( {\frac{{y_{i} }}{{L_{y} }} - \frac{1}{2}} \right)} & {\frac{{y_{i} }}{{L_{y} }} < \frac{1}{2}} \\ {16v_{{{\text{max}}}} \left( {\frac{{y_{i} }}{{L_{y} }} - \frac{1}{2}} \right)\left( {\frac{{y_{i} }}{{L_{y} }} - 1} \right)} & {\frac{{y_{i} }}{{L_{y} }} \ge \frac{1}{2}} \\ \end{array} } \right.. $$

(26)

The maximum velocity along the $x$-axis, represented by $v_{{{\text{max}}}}$, is given by

$$ v_{{{\text{max}}}} = \frac{{L_{y}^{2} F}}{8\mu }. $$

(27)

The system size was set to $L_{x} = 12.2l$ and $L_{y} = 21.1l$, where $l$ is the unit length. Friction from the environment was set to $\xi = 0$. All parameters were non-dimensionalized by units of length, $l$; time, $\tau = 100{\Delta }t$; and energy, $\mu l^{3} /\tau$. The control parameters, $F$, $l_{{\text{d}}}$, and $\lambda$, were explored.

Numerical simulations were performed with several values of $F$, and vertex velocities were observed (Fig. 8b). The resulting velocity profile agreed with the theoretical solution obtained by Eq. (26). By fitting Eq. (26) to the vertex velocities using the least-squares method, the maximum velocity, $v_{{{\text{max}}}}$, was estimated (Fig. 8c). The velocities estimated in the simulations agreed with the theoretical solution. Moreover, the numerical parameters, $l_{{\text{d}}}$ and $\lambda$, were also varied (Fig. 8d, e). Deviations of the numerical solutions from the theoretical solutions given as Eq. (27) were very small and nearly independent of $l_{{\text{d}}}$ and $\lambda$. These results indicate that the planar membrane flow was reproduced in a physically consistent manner. Note that the membrane dynamics are given deterministically in this simulation, whereas the observed velocities were subject to some variability, expressed as standard deviation. This variability was caused by the heterogeneous structure of the triangular mesh, as well as the variation of velocities in the x-axis direction; because the triangular mesh discretizing the membrane is distorted, its structure is heterogeneous, inducing artificial variations in velocities. Moreover, because friction from the environment was set to zero, membrane velocity was determined as a relative velocity, allowing it to physically vary in the x-axis direction.

1.2 A2: Negligible gaps of geometry and energy in dynamic remeshing

Although planar membrane flow was successfully reproduced in the proposed framework (Sect. A1 in Appendix), the flow in 3D space may be affected by the dynamic remeshing algorithm. To clarify whether the artificial effects introduced by the remeshing algorithm are negligible, numerical simulations of rotational cell dynamics were performed under simple shear flow, and gaps in the membrane geometry and energy induced by the remeshing algorithm were evaluated.

The effective energy, $U$ in Eq. (7), is expressed as

$$ U \equiv - nk_{B} T\log \left( {\frac{V}{{V_{{{\text{ref}}}} }}} \right) + \mathop \sum \limits_{i}^{{{\text{vertex}}}} \Gamma_{0} A_{i} + \mathop \sum \limits_{i}^{{{\text{vertex}}}} 2K\frac{{M_{i}^{2} }}{{A_{i} }} $$

(28)

The first term indicates the osmotic energy of the cytoplasm [30], where the constant $nk_{B} T$ is the modulus of the cell volume. Variable $V$ and constant $V_{{{\text{ref}}}}$ are the current and reference volumes of the cell, respectively. The second term in Eq. (28) indicates the surface energy of the membrane, where constant $\Gamma_{0}$ is the surface tension and variable $A_{i}$ is the local surface area around the ith vertex, given by $A_{i} = \mathop \sum \nolimits_{j\left( i \right)}^{{{\text{triangle}}}} a_{j} /3$, where $a_{j}$ is the area of the jth triangle that comprises the ith vertex. Assuming that the first and second terms create a force balance on a spherical geometry when $V = V_{{{\text{ref}}}}$, $nk_{B} T$ may be rewritten as a function of $\Gamma_{0}$ [30]:

$$ nk_{B} T = \Gamma_{0} \left( {\frac{{32\pi V_{{{\text{ref}}}}^{2} }}{3}} \right)^{\frac{1}{3}} $$

(29)

The third term in Eq. (28) indicates the bending energy of the membrane [39], where $K$ is the bending rigidity. Variable $M_{i}$ is the total mean curvature around the ith vertex, given by $M_{i} = \mathop \sum \nolimits_{j\left( i \right)}^{{{\text{edge}}}} l_{j} \theta_{j} /4$, where $l_{j}$ and $\theta_{j}$ are the length and angle of the jth edge connected to the ith vertex, respectively.

A system box with a sufficiently large space was considered in $xyz$-coordinates. To provide a Couette flow as the environment, a simple shear flow was introduced as ${\varvec{u}}_{{\text{e}}} \equiv \left( {qz, 0, 0} \right)$, where $q$ is the gradient in velocity along the $z$-axis. The shear flow exerts shear stress on the cell via the friction force term in Eq. (7), balancing the surface tension in Eq. (28). The physical parameters were non-dimensionalized by length, $l = \left( {3V_{{{\text{ref}}}} /4000\pi } \right)^{1/3}$; time, $\tau = \mu l/\Gamma_{0}$; and energy, $\Gamma_{0} l^{2}$. The remaining parameters were set to $\xi = \mu /2l^{2}$, $K = \Gamma_{0} l^{2}$, $q = 1/\tau$, and ${\Delta }t = 0.01\tau$, while the control parameters, $l_{{\text{d}}}$ and $\lambda$, were explored.

Under shear flow, the cell rotated while maintaining its shape (Fig. 

Fig. 9

Geometric and energetic gaps via remeshing. a. Snapshots of a rotating cell, obtained at $t = 100\tau$. The colors indicate local membrane velocities, ${\varvec{u}}_{i}$. b–e. Average gap ratios of area, volume, curvature, and energy before and after the remeshing operation, as functions of the discretization length $l_{{\text{d}}}$. Results in b–e are from numerical simulations with $\lambda = 0.01$. f–i. Average gaps of area, volume, curvature, and energy before and after the remeshing operation, as functions of the threshold ratio $\lambda$. Results in f–i were obtained from simulations with $l_{{\text{d}}} = l$. Error bars in b–i indicate standard deviations

9a). To evaluate the geometric and energetic gaps, several physical quantities were observed: the volume, $V$; total area, A_t $= \mathop \sum \nolimits_{i} A_{i}$; total mean curvature, $M_{{\text{t}}} = \mathop \sum \nolimits_{i} M_{i}$; and total energy, $U$. During cell rotations, the gaps of $V$, $A_{{\text{t}}}$, $M_{{\text{t}}}$, and $U$ were maintained at low values within certain variances (Fig. 9b–e). Importantly, both the averages and variances of these parameters were reduced by decreasing the discretized length, $l_{{\text{d}}}$ (Fig. 9b–e). Additionally, to evaluate the effects of the correction term in Eq. (13), the cases were compared with and without the correction (Fig. 9b–e). Gaps in the case with the correction were approximately half those without the correction, indicating that the correction term was effective. Moreover, these gaps were reduced by decreasing the threshold ratio, $\lambda$ (Fig. 9f–i). For example, the gaps became less than approximately $0.001\%$ when $l_{{\text{d}}} \le l$ and $\lambda \le 0.01$. These results indicate that the gaps in membrane geometry and energy generated by the remeshing algorithm can be made negligible by modulating $l_{{\text{d}}}$.