Unfitted finite element modelling of surface-bulk viscous flows in animal cells

TL;DR

This work addresses cortex–cytoplasm morphogenesis by formulating a coupled surface-bulk viscous flow model in a fixed Cartesian grid, enabling large deformations without remeshing. It combines a sharp-interface trace FEM for the surface with an aggregated FEM for the bulk, linked through a mechanochemical surface regulator that modulates active tension via diffusion and exchange with the bulk, all tracked by a level-set representation. The main contributions are (i) a stable partitioned unfitted FE framework for surface-bulk pde with robust handling of cut cells, (ii) a time-discrete scheme that couples surface transport, surface and bulk flows, and surface evolution, and (iii) demonstration through axisymmetric and 3D simulations of self-organization, curvature-driven relaxation, and cytokinesis-like deformations. The approach offers a powerful, extendable tool for exploring complex morphogenetic processes in animal cells with high fidelity to moving interfaces and multiscale mechanochemical coupling.

Abstract

This work presents a novel unfitted finite element framework to simulate coupled surface-bulk problems in time-dependent domains, focusing on fluid-fluid interactions in animal cells between the actomyosin cortex and the cytoplasm. The cortex, a thin layer beneath the plasma membrane, provides structural integrity and drives shape changes by generating surface contractile forces akin to tension. Cortical contractions generate Marangoni-like surface flows and induce intracellular cytoplasmic flows that are essential for processes such as cell division, migration, and polarization, particularly in large animal cells. Despite its importance, the spatiotemporal regulation of cortex-cytoplasm interactions remains poorly understood and computational modelling can be very challenging because surface-bulk dynamics often lead to large cell deformations. To address these challenges, we propose a sharp-interface framework that uniquely combines the trace finite element method for surface flows with the aggregated finite element method for bulk flows. This approach enables accurate and stable simulations on fixed Cartesian grids without remeshing. The model also incorporates mechanochemical feedback through the surface transport of a molecular regulator of active tension. We solve the resulting mixed-dimensional system on a fixed Cartesian grid using a level-set-based method to track the evolving surface. Numerical experiments validate the accuracy and stability of the method, capturing phenomena such as self-organised pattern formation, curvature-driven relaxation, and cell cleavage. This novel framework offers a powerful and extendable tool for investigating increasingly complex morphogenetic processes in animal cells.

Figures (13)

• Figure 3: Illustration of an unfitted fe approximation for a bulk pde in $\Omega$. (a) The fe discretisation is built on the $\Omega$-active mesh$\mathcal{T}_h^\Omega$ formed by the grid cells intersecting $\Omega$. The simplest scenario is to approximate the problem with a linear Lagrangian fe space $\mathcal{V}_h^\Omega$, where the shape functions and degrees of freedom (in yellow circles) are associated to the vertices of $\mathcal{T}_h^\Omega$. (b) The fe solution $u_h$ of the problem is found in $\mathcal{T}_h^\Omega$ and its restriction $u_h{_|{_\Omega}}$ is the approximation to the continuous problem. (c) The small cut cell problem. Illustration of a badly cut cell $T$ with a very small cut portion in $\Omega$. If $a$ is a dof of $\mathcal{V}_h^\Omega$, we see that its support inside $\Omega$ is so small that it can lead to ill-conditioning (see Equation \ref{['eq:condition-number']}). (d) A way to fix the small cut cell problem, proposed in, e.g., the Aggregated fem badia2018aggregated, is to extrapolate the dof value $u_a$ in terms of the dof values of an interior cell $T(a)$. Purple arrows indicate the path mapping $a$ to $T(a)$ according to a cell aggregation scheme.
• Figure 4: Schematic representation of the coupled iterative scheme described in Algorithm \ref{['alg:main-loop']}.
• Figure 5: Verification examples.2D axisymmetric surface-bulk fluid sphere problem adapted from happel2012low:(a) Velocity magnitude and streamlines of the steady-state solution. (d) The rate of decay of the approximation error under uniform refinements obeys theoretical estimates Badia2018Mixedjankuhn2021error. 2D axisymmetric dynamic surface:(b) Initial species concentration and velocity fields. Species concentration at the left sextant is 10% higher than elsewhere, setting the system out of equilibrium. (c) Equilibrium shape and solution in the comoving frame at $T = 3.0$. (e) Decay of the experimental errors with uniform refinements. We observe rates of convergence close to quadratic order for all quantities. Note that the results are reflected across the axis of symmetry for illustration purposes.
• Figure 6: Self-organised shape emergence.Tuning the mesh size for analysis of pattern formation. Before the analysis reported in Figure \ref{['fig:pattern']}, we found the minimum mesh size needed to correctly capture the stable to unstable transition in the linear regime, as predicted by Equation \ref{['eq:critical-pe']}. (a) Setting $h = 0.04$ provides a good balance between computational cost and correctly capturing the theoretically predicted appearance of the first instability mode. (b) At $\hbox{Pe} = 11$ and $h = 0.04$, no instabilities appear and concentration heterogeneities slowly diffuse out. (c) At $\hbox{Pe} = 13$ and $h = 0.04$, we observe a paradigmatic first order polar instability forming.
• Figure 7: Self-organised shape emergence.Pattern formation at large hydrodynamic lengths. A perturbation of an homogeneous concentration field leads to spontaneous pattern formation at large enough $\hbox{Pe}$ numbers. (a) Peak concentration and (b) migration velocities for $\hbox{Pe} \in [10,150]$ and $L_\eta \in \{10^3,10^4,10^5\}$. The evolution of the leading mode is coloured along the plotted lines. At large hydrodynamic lengths, we do not observe differences in shapes, flow fields and magnitudes. For this reason, we only plot the results for $L_\eta = 10^4$. (c) Steady-state results obtained for leading modes 1-2-3 at, respectively, $\hbox{Pe} = 30,60,120$.