††thanks: These authors contributed equally.††thanks: These authors
contributed equally.

# Dimensionality reduction in bulk–boundary reaction–diffusion systems

Tom Burkart  Arnold Sommerfeld Center for Theoretical Physics and Center for
NanoScience, Department of Physics, Ludwig-Maximilians-Universität München,
Theresienstraße 37, D-80333 München, Germany  Benedikt J. Müller  Arnold
Sommerfeld Center for Theoretical Physics and Center for NanoScience,
Department of Physics, Ludwig-Maximilians-Universität München, Theresienstraße
37, D-80333 München, Germany  Erwin Frey<EMAIL_ADDRESS>Arnold Sommerfeld Center
for Theoretical Physics and Center for NanoScience, Department of Physics,
Ludwig-Maximilians-Universität München, Theresienstraße 37, D-80333 München,
Germany  Max Planck School Matter to Life, Hofgartenstraße 8, D-80539 München,
Germany

(August 28, 2024)

###### Abstract

Intracellular protein patterns regulate many vital cellular functions, such as
the processing of spatiotemporal information or the control of shape
deformations. To do so, pattern-forming systems can be sensitive to the cell
geometry by means of coupling the protein dynamics on the cell membrane to
dynamics in the cytosol. Recent studies demonstrated that modeling the
cytosolic dynamics in terms of an averaged protein pool disregards possibly
crucial aspects of the pattern formation, most importantly concentration
gradients normal to the membrane. At the same time, the coupling of two
domains (surface and volume) with different dimensions renders many standard
tools for the numerical analysis of self-organizing systems inefficient. Here,
we present a generic framework for projecting the cytosolic dynamics onto the
lower-dimensional surface that respects the influence of cytosolic
concentration gradients in static and evolving geometries. This method uses a
priori physical information about the system to approximate the cytosolic
dynamics by a small number of dominant characteristic concentration profiles
(basis), akin to basis transformations of finite element methods. As a proof
of concept, we apply our framework to a toy model for volume-dependent
interrupted coarsening, evaluate the accuracy of the results for various basis
choices, and discuss the optimal basis choice for biologically relevant
systems. Our analysis presents an efficient yet accurate method for analysing
pattern formation with surface–volume coupling in evolving geometries.

##  I Introduction

Living organisms ensure their viability by precisely orchestrating a wide
array of cellular functions, ranging from subcellular information processing
to morphogenesis. At the heart of these functions lie out-of-equilibrium
molecular systems: proteins whose spatio-temporal organization in cells is
driven by chemical interactions with other proteins and transport across the
cell, commonly referred to as reaction–diffusion systems. From these two
building blocks, complex information processing systems can emerge that allow
the transmission of chemical signals within the cell or encode temporal and
spatial cues [1, 2, 3]. Popular examples of such pattern-forming protein
systems are the Min system of Escherichia coli which exhibits pole-to-pole
oscillations in vivo [4] and spiral waves or labyrinth patterns in vitro [5];
localization of the budding site via Cdc42 in Saccharomyces cerevisiae [6]; or
polarity establishment by PAR proteins in Caenorhabditis elegans [7].

Figure 1: Schematic visualization of the projection of a time-dependent bulk
field $c(\bm{s},z,t)$ (indicated by green (gray) slice) in a three-dimensional
bulk–boundary system onto the two-dimensional boundary of the domain via an
approximation using multiple characteristic profiles $v_{k}(z)$ and
corresponding time dependent membrane fields $u_{k}(\bm{s},t)$ (indicated by
green (gray) lines).

In dynamically deforming geometries, pattern-forming systems may gain the
ability to engage in mechanochemical feedback loops that can give rise to even
more interesting self-organizing dynamics and that play a key role in the
regulation of many cellular functions. Examples on the sub-cellular scale
include the sensing and generation of membrane curvature by membrane-binding
proteins [8, 9, 10, 11, 12] or the adaptive establishment of cell polarity
[13]. On the single-cell and tissue level, pattern formation in deforming
geometries has been appreciated in the context of biochemical coordination of
cell motility [14, 15, 16], cell shape changes [17, 18, 19, 20], and the
control of tissue morphogenesis [21, 22]. A central common feature of many
pattern-forming systems in dynamic geometries is the coupling of two
concurrently deforming domains – the cell volume (or bulk) and the cell
membrane (boundary) – which models need to consider explicitly in order to
accurately capture the mechanochemical coupling of the protein dynamics to the
geometry [3].

The realization of such systems in numerical simulations, however, offers two
fundamental challenges: first, most numerical frameworks are not designed for
solving reaction–diffusion dynamics in a deforming geometry, and including
such deformations often abates the performance considerably. Second, since the
diffusion of membrane-bound proteins is typically slow compared to cytosolic
proteins [23], the length scale of protein patterns on membranes can be by
orders of magnitude shorter than the length scale of protein gradients in the
cell volume or the length scale of deformations. This requires a high spatial
resolution of the membrane and its vicinity, resulting in a fine mesh with
high number of degrees of freedom that further slows down simulations.

A plethora of mathematical studies have focused on reaction–diffusion dynamics
on the surface of a deforming geometry. Popular methods that are used in such
cases include specialized applications of the finite element method (FEM) [24,
25, 26], level-set methods [27, 28], and mesh-free approaches [29, 30].
However, it is often challenging to generalize these approaches to systems
with bulk–boundary coupling, i.e., systems with dynamics both in the volume
and on the surface of the deforming geometry, as this requires handling the
deformation effects in two different domains simulatneously and in a self-
consistent manner [31]. Many studies of biological systems handle this
limitation by averaging the volume dynamics either over infinitesimal surface
patches (projection) or over the entire volume (reservoir) [32, 33]. While
this approach yields acceptable results when the volume dynamics are purely
diffusive [34], it was recently demonstrated that accounting for the bulk
dynamics is crucial for accurate predictions in the presence of bulk reactions
[35, 36, 37]. This is because bulk reactions generically lead to protein
concentration gradients in the cytosol which, by virtue of the bulk–boundary
coupling, directly affect the protein dynamics on the membrane and thereby
provide the basis for one type of geometry sensing [3]. In studies committed
to faithfully representing the reaction–diffusion dynamics in deforming
geometries, the phase–field method has emerged as a promising strategy [38,
39, 40, 41, 42, 43, 44, 15, 45, 46]. Phase–field models represent the dynamic
geometry as an indicator function, allowing for arbitrary shape changes and
even topological changes [47, 39]. However, this benefit comes at the cost of
either requiring fine mesh resolution everywhere in the simulated domain or
using adaptive mesh refinement [48, 49]. Furthermore, coupling bulk dynamics
to surface dynamics in phase–field models raises additional challenges that
are subject of ongoing development [50, 51].

Here, we present an approach that exploits a priori knowledge about cytosolic
gradients for modeling reaction–diffusion dynamics with bulk–boundary coupling
in deforming geometries. Rather than solving the ensuing partial differential
equations (PDEs) in a meshed bulk, we propose to project the bulk dynamics
onto the surface of the deforming geometry by decomposing the bulk dynamics
into a few dominant basis functions, essentially reducing the spatial
dimension of the system by one [Fig. 1]. We show how to design projection
methods based on information about the steady state of the system at hand and
compare three different techniques – based on step functions, polynomials, and
exponential functions – with respect to their accuracy in flat and in
deforming geometries. Our analysis demonstrates that good approximations of
the actual dynamics can be achieved already with a single nonlinear basis
function, whereas a naive averaging of the bulk dynamics often fails to
capture the relevant dynamics entirely. Furthermore, we find that the
projection method outperforms standard FEM approaches by an up to five-fold
speedup in computation time, with most significant improvements in systems
with large volume.

This paper is structured as follows: In Sec. II, we first introduce a generic
projection method in a static one-dimensional geometry and subsequently
generalize our approach to dynamic higher-dimensional geometries. We show how
to transfer techniques from FEM implementation to arbitrary basis choices and
how a PDE defined on a volume domain can be approximated by a set of PDEs
defined only on the volume’s boundary. In Sec. III we test our approach on a
toy model that shows geometry-dependent arrested coarsening and compare
different projection methods based on the quality of the approximation and the
computational efficiency. We conclude with a concise summary and an outlook.

##  II Bulk projection

In this section, we introduce a method to project the dynamics of proteins in
the bulk onto the membrane for a flat geometry. We illustrate this method
using a generic model for a single protein species that can bind to and unbind
from the membrane and are degraded linearly in the bulk. We then extend the
projection method to dynamically deforming membranes.

###  II.1 Static geometry

Figure 2: Schematic visualization of flat geometries. (a) A one-dimensional
bulk domain $\mathcal{B}$ (cytosol, blue/gray) of size $h$ is bounded by a
reactive surface $\mathcal{S}$ (membrane, filled circle) on one side and by a
no-flux boundary on the other side (open circle). (b) For a two-dimensional
bulk domain of height $h$ and width $L$, the surface $\mathcal{S}$ is a one-
dimensional line.

Consider a one-dimensional bulk $\mathcal{B}$ of height $h$ [Fig. 2a]. The
bulk (cytosolic) protein concentration on this line, denoted by $c(z,t)$, is
assumed to obey no-flux (Neumann) boundary conditions at the top of the line
($z=h$) and Robin boundary conditions – where reactive fluxes are balanced by
diffusive fluxes onto the boundary – at the bottom of the line (membrane at
$z=0$, denoted by $\mathcal{S}$). The protein concentration on the membrane is
denoted by $m(t)$.

A generic (non–mass-conserving) bulk–boundary reaction–diffusion system for a
single protein species with linear degradation in the bulk can be written as

$\displaystyle\partial_{t}c(z,t)$
$\displaystyle=D_{c}\,\partial_{z}^{2}c(z,t)-\lambda\cdot c(z,t)\,,$ (1a)
$\displaystyle\partial_{t}m(t)$ $\displaystyle=f_{0}(m(t),c(0,t))\,,$ (1b)

where the bulk dynamics are coupled to the membrane via Robin and no-flux
boundary conditions

$\displaystyle-D_{c}\,\left.\partial_{z}c(z,t)\right|_{z=0}$
$\displaystyle=-f_{0}(m(t),c(0,t))\,,$ (2a) $\displaystyle
D_{c}\,\left.\partial_{z}c(z,t)\right|_{z=h}$ $\displaystyle=0\,.$ (2b)

Here, $D_{c}$ denotes the diffusion constant in the bulk, $\lambda$ is a
cytosolic degradation rate, and $f_{0}(m(t),c(0,t))$ denotes the reactive
fluxes at the surface $\mathcal{S}$. A generic choice for the reaction at the
boundary is

$f_{0}(m(t),c(0,t))=a(m(t))\cdot c(0,t)-d(m(t))\cdot m(t)\,,$

where the attachment $a$ and detachment $d$ to and from the surface
$\mathcal{S}$ may include (autocatalytic) nonlinear interactions.

To project the dynamics in a one-dimensional bulk onto the zero-dimensional
surface $\mathcal{S}$, we aim to rewrite the partial differential equation
(PDE) for the bulk dynamics in Eq. (1) as a set of $N$ conveniently chosen
ordinary differential equations (ODEs) that approximate the exact solution.
Analogous to reaction–diffusion equations that can be split into reaction and
diffusion parts [Eq. (1)] and boundary conditions [Eq. (2)], we thus aim for a
set of ODEs in the form of

$\partial_{t}u_{k}(t)=g_{\text{React}}+g_{\text{Diff}}+g_{\text{BC}}$ (3)

for $k\in\\{0,\ldots,N{-}1\\}$. The symbolic functions
$g=g\bigl{(}\\{u_{k}(t)\\},m(t)\bigr{)}$ are placeholders for the
contributions to the ODE stemming from the bulk reactions, bulk diffusion, and
the boundary conditions, respectively. How can one derive such a set of ODEs
(or, in higher dimensions, PDEs)?

To answer this question, we draw inspiration from numerical methods for
solving PDEs, specifically the _finite element method_ (FEM). The general idea
of FEM approaches is to approximate the exact dynamics by a linear combination
of contributions derived for simple basis functions (finite elements), where
each basis function represents the field value at a specific point in the
meshed simulation domain [52, 53]. Since, for generic applications, the
dynamics of the concentration fields are not known a priori, it is vital for
an efficient FEM solver to use simple basis functions with little overlap on a
well-resolved mesh. For physical systems, however, certain properties of the
field dynamics can be derived a priori. For example, in a reaction–diffusion
system with linear bulk reactions as in Eq. (1), the characteristic diffusive
length scale in the bulk is $\ell=\sqrt{D_{c}/\lambda}$ [54]. In the
following, we show how a more sophisticated choice of basis functions that
exploits such knowledge can simplify the problem and reduce the complexity of
numerical implementations solving the system’s dynamics.

Figure 3: Comparison of different basis choices. Individual basis functions
are indicated by different shades of gray. (a) Typical choice of piecewise
linear basis functions used in FEM with irregular mesh (mesh points indicated
by open circles). (b) Orthonormalized basis derived from
$v_{k}=\cosh\bigl{(}(h{-}z)/\ell_{k}\bigr{)}$ with
$\ell_{k}=\\{\infty,1,2\\}$. This basis choice can capture strong gradients
close to the membrane ($z=0$) paired with weak gradients at large heights, but
not vice versa. (c) Orthonormalized basis derivd from $v_{k}=(z-h)^{2k}$. (d)
Orthonormalized piecewise constant basis with increasing layer heights $h_{k}$
for fine resolution of the membrane-proximal region.

While FEM in the context of numerical simulations typically uses a single
basis function per mesh point in the simulated domain [53], we pursue a
different approach using multiple basis functions $v_{k}(z)$ anchored to the
boundary of the domain [Fig. 3]. The different basis functions are designed to
reflect the geometry of the bulk domain and to ensure that the relevant
aspects of concentration gradients perpendicular to the boundary can be
captured appropriately. In contrast to FEM, the coefficients $u_{k}(t)$ are
defined for mesh points on the membrane so that all bulk variations are
captured by the height-dependent basis functions $v_{k}(z)$. With this
(incomplete) basis, the actual bulk dynamics are approximated by

$c(z,t)\approx\sum\limits_{k}u_{k}(t)\cdot v_{k}(z)\,.$ (4)

Inserting this ansatz into the reaction–diffusion equations (1) and (2)
yields, after partially integrating the diffusion term and taking the inner
product $\langle\cdot,\cdot\rangle$ with a basis function $v_{k}$ (Galerkin
projection [53]), a set of coupled ODEs for the coefficients $u_{k}$ (detailed
derivation in Appendix A):

$\displaystyle G_{kl}\,\partial_{t}u_{k}(t)$
$\displaystyle=\underbrace{-f_{0}\bigl{(}m(t),\\{u_{k}(t)\\}\bigr{)}\,v_{l}(0)}_{g_{\text{BC}}}+$
$\displaystyle\qquad\underbrace{-D_{c}\,A_{ml}\,u_{m}(t)}_{g_{\text{Diff}}}-\underbrace{\lambda\,G_{lm}\,u_{m}(t)}_{g_{\text{React}}}\,.$
(5)

Here and in the following, we use Einstein sum convention over double indices.
With the inner product

$\langle v(z),w(z)\rangle=\int_{0}^{h}\\!\\!\\!\mathrm{d}z\;v(z)\,w(z)$

the _mass matrix_ (or Gram matrix) $G_{kl}$ and _stiffness matrix_ $A_{kl}$
are defined as

$\displaystyle G_{kl}$ $\displaystyle=\langle v_{k}(z),v_{l}(z)\rangle\,,$
$\displaystyle A_{kl}$
$\displaystyle=\langle\partial_{z}v_{k}(z),\partial_{z}v_{l}(z)\rangle\,.$ (6)

From here on we will omit all function arguments unless needed for clarity.

The set of ODEs in Eq. (5), together with the basis $\\{v_{k}\\}$, yield an
approximation of the time evolution of the bulk field $c$, where the quality
of the approximation critically depends on the basis choice. What is a good
choice for the basis $\\{v_{k}\\}$? In general, the basis needs to be chosen
such that it can capture the main features of the bulk dynamics. Specifically,
this requires that both the field $c$ as well as its gradient $\partial_{z}c$
are approximated well by the projection $u_{k}\,v_{k}$. Indeed, it has been
shown that the coefficients $u_{k}$ obtained from $u_{k}=(G^{-1})_{kl}\langle
v_{l},c\rangle$ yield the best possible approximation for the field $c$ for a
given basis choice $\\{v_{k}\\}$ [53]. However, no such statement exists for
the gradient $\partial_{z}c$, and the same coefficients $u_{k}$ can result in
a low-quality approximation of the gradient for improper basis choices. Here
we present two ways to resolve this problem:  (i) by choosing a basis that is
expected a priori to yield good approximations for both field $c$ and gradient
$\partial_{z}c$, or  (ii) by making additional choices for the gradients and
thereby correct for the low-quality approximation of the gradients.

The former approach strongly depends on the physical problem that is being
studied. For example, for protein reaction–diffusion systems with linear bulk
reactions as specified in Eq. (1) and no-flux boundary conditions at $z{=}h$,
the steady state distribution can be derived as
$c(z)\sim\cosh\bigl{(}(h-z)/\ell\bigr{)}$, with $\ell=\sqrt{D_{c}/\lambda}$
[54]. In this case, a natural choice for the basis would be a set of
hyperbolic cosine functions on varying length scales [Fig. 3b], e.g.,

$\\{v_{k}(z)\\}=\left\\{1,\,\cosh\left(\frac{h-z}{\ell}\right),\,\cosh\left(\frac{h-z}{2\ell}\right),\,\ldots\right\\}\,.$

A mathematically more tractable choice for the same system is
$\\{v_{k}(z)\\}=\\{(h{-}z)^{k}\\}$, which corresponds to a polynomial fit to
the field $c$ and can be expected to approximate shallow gradients well [Fig.
3c].

For the latter approach revolving around making additional choices for the
gradients, we only discuss a basis composed of orthogonal step functions, as
this choice will prove highly useful for deforming geometries later on. For
this, the bulk is divided in $N$ sections of height $h_{k}$, where each basis
function $v_{k}$ takes the form of a rectangle function [Fig. 3d]. Since the
derivative of these basis functions at the interfaces $z_{k}$ are not well-
defined, we instead approximate the gradients in between the centers of the
step functions as an additional set of piecewise constant functions, as
derived in Appendix B. To highlight this additional approximation, we indicate
auxiliary stiffness matrices constructed from such an artificial gradient
choice by a bar over the symbol ($\bar{A}_{kl}$). For piecewise constant
gradients, the resulting set of ODEs for the coefficients $u_{k}$ then reads
[Appendix B]

$\displaystyle\partial_{t}u_{k}$
$\displaystyle=-\delta_{k0}\,f_{0}\,v_{0}(0)-D_{c}\,\bar{A}_{kl}\,u_{l}+\lambda\,u_{k}\,,$
(7a) $\displaystyle\bar{A}_{kk}$
$\displaystyle=\frac{2}{h_{k}+h_{k{+}1}}+\frac{2}{h_{k}+h_{k{-}1}}\,,$ (7b)
$\displaystyle\bar{A}_{k,k{\pm}1}$
$\displaystyle=-\frac{2}{h_{k}+h_{k{\pm}1}}\,,$ (7c)

where $\delta_{kl}$ is the unit matrix and all other entries $\bar{A}_{kl}=0$.
It is instructive to compare this choice with standard finite element methods
on a predefined mesh: in FEM implementations, basis functions are typically
non-zero only in the direct vicinity of a specific mesh point, similar to the
step functions that are non-zero only for a fraction of the bulk domain [Fig.
3d] [53]. Furthermore, FEM basis functions are typically piecewise linear,
such that their gradients are piecewise constant and thus the stiffness matrix
takes a similar form as in Eq. 7a. Rather than fully committing to the FEM
approach, we here stick to piecewise constant basis functions and the
auxiliary stiffness matrix $\bar{A}_{kl}$ since this will greatly simplify the
generalization to deforming geometries in the following section.

###  II.2 Dynamically deforming geometries

Figure 4: Schematic visualization of a deforming geometry. The bulk domain
$\mathcal{B}$ (cytosol, blue/gray) is parametrized by a vector
$\bm{r}(s,z;t)$, with $z$ the distance from the reactive boundary
$\mathcal{S}$ (membrane, thick black line) measured along the normal vector
$\bm{n}(s,t)$. The boundary itself is parametrized by
$\bm{r}_{\mathcal{S}}(s,t)=\bm{r}(s,0;t)$.

We now generalize the idea of bulk projection to dynamically deforming
geometries. Consider a two-dimensional domain $\mathcal{B}$ parametrized by
$\bm{r}(s,z)$, where $s$ denotes the (arc-length) parametrization of the one-
dimensional membrane $\mathcal{S}$, $\bm{r}_{\mathcal{S}}(s,t)$, and
$z\in[0,h]$ is the distance from the membrane [Fig. 4]:

$\bm{r}(s,z;t)=\bm{r}_{\mathcal{S}}(s,t)+z\,\bm{\hat{n}}(s,t)\,.$ (8)

The vector $\bm{\hat{n}}$ is the normal vector on the membrane. The bulk
dynamics specified in Eq. (1) in such a dynamic geometry read [12]

$\displaystyle\frac{1}{\sqrt{g}}\partial_{t}\bigl{(}\sqrt{g}\,c\bigr{)}$
$\displaystyle=D_{c}\,\Delta_{\text{LB}}c+\lambda\cdot c\,,$ (9a)
$\displaystyle\Delta_{\text{LB}}c$
$\displaystyle=\frac{1}{\sqrt{g}}\partial_{i}\bigl{(}\sqrt{g}g^{ij}\partial_{j}c\bigr{)}\,,$
(9b)

with the Laplace-Beltrami operator $\Delta_{\text{LB}}$, the metric tensor
$g_{ij}$, and the square root of the metric determinant
$\sqrt{g}\equiv\sqrt{\det(g)}$. Furthermore, we generalize the inner product,
the Gram matrix $G$, and the stiffness matrix $A$ to the deformable geometry:

$\displaystyle\langle u(z),v(z)\rangle^{ij}$
$\displaystyle=\int_{0}^{h}\\!\\!\\!\mathrm{d}z\;\sqrt{g(s,z;t)}\,g^{ij}(s,z;t)\,u(z)\cdot
v(z)\,,$ $\displaystyle G_{kl}^{ij}(s,t)$
$\displaystyle=\langle\,v_{k},v_{l}\rangle^{ij}\,,$ (10) $\displaystyle
A_{kl}^{ij}(s,t)$
$\displaystyle=\langle\partial_{z}v_{k},\partial_{z}v\rangle^{ij}\,,$

so that $G_{kl}\equiv G^{zz}_{kl}=\langle v_{k},v_{l}\rangle^{zz}$, and
analogous for $A_{kl}$. In Appendix A, we use these definitions to derive the
weak formulation and Galerkin projection [53] of the bulk dynamics in a
deforming geometry as specified in Eq. (8). The resulting generalized set of
PDEs for the coefficients reads

$\displaystyle G_{kl}\partial_{t}u_{k}$
$\displaystyle=\underbrace{-u_{m}\,(\partial_{t}G)_{lm}}_{g_{\text{Geom}}}+\underbrace{D_{c}\,\partial_{s}\bigl{(}G^{ss}_{lm}\,\partial_{s}u_{m}\bigr{)}}_{g_{\mathcal{S}}}+$
(11)
$\displaystyle\quad+D_{c}\bigl{[}\sqrt{g}\partial_{z}c\,v_{l}\bigr{]}_{0}^{h}-D_{c}\,u_{m}\,A_{lm}+\lambda\,u_{m}\,G_{lm}\,.$

In comparison to the set of ODEs in Eq. (5), two new terms have emerged in
addition to the redefined coupling matrices: The contribution
$g_{\text{Geom}}$ accounts for deformations of the domain (dilation and
compression), and $g_{\mathcal{S}}$ accounts for the diffusion parallel to the
membrane.

Strikingly, this formulation holds for any basis where all elements
$\\{v_{k}\\}$ are weakly differentiable on $[0,h]$. In addition, the
formulation can be adjusted for a basis of orthogonal step functions by
appropriate rescaling of the stiffness matrix $\bar{A}_{kl}$ [Appendix B].
This offers significant advantages compared to an explicit simulation of the
bulk dynamics: first, the number of degrees of freedom can be reduced
significantly by avoiding a meshed bulk domain. Second, similar to other
methods that use an explicit parametrization of the geometry, our approach
also makes a dynamic remeshing of the deforming bulk obsolete. These
advantages come at the cost of having to settle on a basis a priori, where the
quality of the approximation critically depends on the basis choice.

With the projection step in Eq. (11), one therefore converts a
reaction–diffusion system with coupled PDEs in the bulk and on the membrane
into a set of PDEs that are defined only on the membrane. This approach yields
an effective reduction of the system’s spatial dimension (e.g., from a three-
dimensional bulk to a two-dimensional membrane) while accounting for bulk
gradients in an approximate manner and thereby respects the effect of the bulk
geometry on the protein dynamics. The accuracy of this approximation depends
on the basis choice $\\{v_{k}\\}$. These benefits come at the cost of
introducing one additional PDE on the membrane for each basis function.

##  III Examples

Figure 5: (a) Schematic reactions in the geometry-sensitive coarsening toy
model. A protein can bind to an unbind from the membrane in its active state
(solid). In the bulk, the protein undergoes deactivation (dashed). On the
membrane, inactive proteins can be activated assisted by active bulk proteins.
The inactive proteins are assumed to be abundant (reservoir) and undergo
linear binding to and unbinding from the membrane, thereby maintaining a
constant concentration on the membrane. (b) The steady-state concentration of
the bulk field close to the membrane $c^{*}_{\mathcal{S}}$ increases with bulk
height for $h\lesssim\ell$ but is approximately constant for $h\gg l$.

In the following, we present and discuss an explicit application of our
dimensionality reduction method for a simple toy model, with particular focus
on comparing different basis choices and on the influence of the bulk
geometry. For this, we extend a previously studied model showing interrupted
coarsening by introducing an explicit sensitivity to the bulk geometry [55,
56]. This choice allows to easily quantify the effect of the bulk geometry
and, most importantly, how well the bulk dynamics are captured by the
projection method by means of the pattern length scale after coarsening is
interrupted. In this model, we consider a protein that switches between an
active and an inactive state and prevails either in the bulk or on the
membrane [Fig. 5a]. In the active state, the protein can bind and recruit
itself to the membrane, and undergoes enzyme-mediated detachment. Active
proteins deactivate both on the membrane and in the bulk, but are assumed to
be reactivated only on the membrane and assisted by membrane-proximal active
bulk proteins. For simplicity, we furthermore assume that the inactive
proteins comprise an abundant reservoir maintaining a constant concentration
on the membrane and are therefore not modeled explicitly. The
reaction–diffusion equations for the active proteins on the membrane $m(x,t)$
and in the bulk $c(\bm{x},t)$ then read [55, 56]

$\displaystyle\partial_{t}m$
$\displaystyle=D_{m}\,\partial_{x}^{2}m+f(m,\left.c\right|_{\mathcal{S}})+\zeta(m,\left.c\right|_{\mathcal{S}})\,,$
(12a) $\displaystyle\partial_{t}c$
$\displaystyle=D_{c}\,\nabla^{2}c-\lambda\,c\,,$ (12b)

with Robin boundary conditions

$D_{c}\,\left.\bm{n}\cdot\nabla
c\right|_{\mathcal{S}}=-f(m,\left.c\right|_{\mathcal{S}})\,,$ (13)

where the mass-conserving binding kinetics of active proteins are given by

$f(m,\left.c\right|_{\mathcal{S}})=(1+m)\cdot\left.c\right|_{\mathcal{S}}-\frac{m}{1+m}\,.$
(14)

Here, $\left.c\right|_{\mathcal{S}}=c(x,z{=}0)$ denotes the cytosolic
concentration at the membrane, and we have set all reaction rates and
parameters to $1$ for simplicity. In addition, the non-mass-conserving part
corresponding to the (de-)activation of proteins on the membrane (by exchange
with an abundant reservoir [Fig. 5a]) is

$\zeta(m,\left.c\right|_{\mathcal{S}})=p\,\left.c\right|_{\mathcal{S}}-\epsilon\,m$
(15)

with an effective activation rate $p$ and deactivation rate $\epsilon$. With
Robin boundary conditions at the membrane and linear degradation, the bulk
concentration has an associated length scale $\ell=\sqrt{D_{c}/\lambda}$
characterizing gradients at the steady state [54]. For simplicity, we further
nondimensionalize this toy model by enforcing $\ell=1$, so that the system
height $h$ as well as all other distances will always be given in units of the
characteristic length scale (full parameter list in Table 1).

In Fig. 6a, we show kymographs of two realizations of this system on a domain
of width $L=20$ and at bulk heights $h=1$ and $h=10$ obtained from FEM
simulations that explicitly account for the bulk dynamics. The system
undergoes interrupted coarsening, where initially multiple peaks of high
protein concentration on the membrane (green) form that subsequently merge or
vanish (_coarsening_) until only a few peaks remain (_interruption_). After
coarsening is interrupted, peaks rearrange to accomodate approximately
equidistant spacing $\Lambda$ between them. A signature of interrupted
coarsening is the mean peak distance $\mu_{\Lambda}$ saturating at finite
values [Fig. 6b]. For the reaction–diffusion system in Eq. (12), the
coarsening process is predominantly controlled by the (de-)activation of
proteins on the membrane $\zeta(m,\left.c\right|_{\mathcal{S}})$ [56] and thus
depends on the membrane-proximal bulk concentration
$\left.c\right|_{\mathcal{S}}$. Since this membrane-proximal bulk
concentration in turn depends on the system height $h$ [Fig. 5b], the
coarsening interruption is sensitive to the bulk geometry. Importantly,
similar to the steady-state concentration $c^{*}_{\mathcal{S}}(h)$ shown in
Fig. 5b, we expect a nonlinear dependence of the mean peak distance on the
system height $\mu_{\Lambda}(h)$, with an approximately linear dependence for
$h\lesssim\ell$ and no variation for $h\gg\ell$.

Figure 6: (a) Typical kymographs $m(x,t)$ of the model specified in Eqs. (12)
for $h=1$ (left column) and for $h=10$ (right column). The system undergoes
interrupted coarsening and subsequent wavelength selection. Time is shown on a
log scale. The peak distance $\Lambda$ is indicated in the inset. (b) Mean
peak distance $\mu_{\Lambda}$ for the two kymographs shown in panel (a). At
large system heights, coarsening is interrupted at considerably higher mean
peak distances than for small system heights. Figure 7: Approximation of the
steady-state bulk concentration profile $c(z)$ for the system specified in Eq.
(12) (black, dashed) using finite bases $v_{k}(z)$ for bulk height $h=2\ell$.
Top row shows the best approximation by least squares $u(z)=u_{k}\,v_{k}(z)$
(green/gray), bottom row shows the deviation from the exact profile $\delta
u(z)=c(z)-u(z)$. (a) Naive projection averaging over the bulk height ($N=1$).
(b) A basis using a hyperbolic cosine fully captures the steady-state
concentration profile ($N=1$). (c) A polynomial basis yields good
approximations for sufficiently small bulk heights ($N=2$). (d,e) A piecewise
constant basis yields considerable deviations from the exact profile which can
be reduced by increasing the number of basis functions/“layers” ($N\geq 2$).

###  III.1 Basis choices

As a naive direct projection of the bulk dynamics onto the membrane one may
choose to average the bulk field $c$ over the $z$-direction. This corresponds
to choosing a basis $v_{k}(z)=\\{1\\}$ with a constant function as the only
basis element [Fig. 7a]. By definition, this ansatz disregards all gradients
in the bulk field. These gradients are shallow for small bulk heights
$h<\ell$, however they become significant when the bulk height is larger than
the characteristic length scale of bulk gradients, $h>\ell$. Consequently,
this naive projection captures the average pattern length scale
$\mu_{\Lambda}$ only for small bulk heights but fails for large bulk heights
[Fig. 8a].

In the next step, we extend the projection basis by one additional element to
capture the bulk gradients, i.e., variations in $z$-direction. To specify this
additional basis element, we use the system’s laterally homogeneous steady
state (constant parallel to the membrane), which can be calculated
analytically for the linear reactions as defined in Eq. (12) [54]:

$c^{*}(\bm{x})=c^{*}(z{=}0)\cdot\frac{\cosh\bigl{(}\frac{h-z}{\ell}\bigr{)}}{\cosh(h/\ell)}\,,$
(16)

where $c^{*}(z{=}0)=c^{*}_{\mathcal{S}}$ is the steady state concentration at
the membrane. In the rescaled system with $\ell=1$, the steady state profile
can be represented exactly by a one-dimensional basis
$v_{k}(z)=\\{\cosh(h{-}z)\\}$. Using this basis, gradients on the length scale
of the characteristic scale $\ell$ are captured exactly [Fig. 7b].

In between these two limiting cases (naive averaging using $v_{k}=\\{1\\}$ and
full recovery of the steady state using $v_{k}=\\{\cosh(h{-}z)\\}$) a plethora
of other basis choices can be imagined. The suitability of a specific basis
choice depends on the desired application. For example, while a hyperbolic
cosine-derived basis captures the steady state bulk profile exactly in a flat
geometry, no closed-form expressions for required Gram tensor $G^{ij}_{kl}$
and stiffness tensor $A^{ij}_{kl}$ as defined in Eq. (10) may exist for
deformed geometries with non-trivial metric $g_{ij}$, making them unsuitable
for numerical implementations. In particular, it is desirable to calculate the
Gram and stiffness tensors as functions of time – which requires the inner
products to be solved analytically – since this allows a more flexible
implementation compared to recalculating the tensors after each time step by
solving the inner products numerically. To achieve this, one may expand the
exact profile as a power series that respects the no-flux boundary condition
at $z{=}h$, i.e., using a polynomial basis with even powers [Fig. 7c].

Alternatively, one may choose a set of piecewise constant basis functions that
splits the bulk into $N$ distinct sections of height $h_{k}$ to approximate
the bulk concentration profiles [Fig. 7d]. Similar to finite element methods,
the coefficients $u_{k}$ then denote the average value of the bulk field in
the $k^{\text{th}}$ section. By tuning the section heights $h_{k}$ the
approximation can be refined in regions with comparably sharper bulk
gradients. For example, the membrane-proximal region for the toy model in Eq.
(12) requires a finer resolution along the $z$-direction than regions at the
opposing no-flux boundary due to the shallow gradient at $z{=}h$ [Fig. 7d].
The accuracy of the approximation can further be improved by increasing the
number of basis functions [Fig. 7e].

The comparison of the approximated steady-state concentration profile with the
exact counterpart can provide preliminary intuition about appropriate basis
choices. However, when studying the dynamic case, additional aspects of the
bulk concentration may become relevant. For example, the hyperbolic cosine in
Fig. 7b can capture the long-term (steady state) dynamics but will not account
for sharp bulk gradients as the steady state is approached. How do the basis
choices perform for approximating the protein _dynamics_?

###  III.2 Flat and static geometry

To address this question, we first focus on flat and static geometries [Fig.
2]. Since the projection of the bulk dynamics onto the membrane is a
qualitative approximation, it is not expected that patterns obtained from the
full system and from the approximation match quantitatively. Instead, we are
interested in capturing the qualitative characteristics of the patterns and
compare the pattern’s statistical properties to measure the quality of a basis
choice. For the interrupted coarsening system in Eq. (12), the most important
statistical property is the mean value $\mu_{\Lambda}$ of the distance between
peaks $\Lambda$ after coarsening is interrupted [Fig. 6].

For this, we solved the PDEs for the full system and for the projected
counterparts on a domain of length $L=50$ with varying bulk height
$10^{-1}\leq h\leq 10^{1}$ until coarsening was interrupted (total simulation
time $T=10^{7}$; random perturbations around the homogeneous steady state as
initial conditions) and extracted the mean peak distance $\mu_{\Lambda}$ at
the final time point. To reduce artifacts from initial conditions, the results
are averaged over five repeats with varying initial conditions each. In Fig.
8a, we compare these results for the full system (_reference_) as well as
three different basis choices derived from hyperbolic cosines:

$\displaystyle N$ $\displaystyle=1:$ $\displaystyle v_{k}$
$\displaystyle=\bigl{\\{}\cosh(h{-}z)\bigr{\\}}\,,$ $\displaystyle N$
$\displaystyle=2:$ $\displaystyle v_{k}$
$\displaystyle=\bigl{\\{}1,\cosh(h{-}z)\bigr{\\}}\,,$ $\displaystyle N$
$\displaystyle=3:$ $\displaystyle v_{k}$
$\displaystyle=\bigl{\\{}1,\cosh(h{-}z),\cosh\bigl{(}\tfrac{h{-}z}{\ell_{2}}\bigr{)}\bigr{\\}}$

with a variable length scale $\ell_{2}$ in the latter case. All basis choices
reliably capture the characteristic shape of the reference $\mu_{\Lambda}(h)$,
however they consistently overestimate the true mean peak distance [Fig. 8a].
Importantly, even though the basis with $v_{k}^{(N{=}1)}=\\{\cosh(h{-}z)\\}$
can capture the steady state exactly as shown in Fig. 7b, additional basis
functions are required to reduce the deviation from the reference in the
dynamic case.

Figure 8: Comparison of the mean peak distance $\mu_{\Lambda}$ after
coarsening using the projection method with different basis choices. Reference
values obtained from explicit FEM simulations including the bulk are shown as
dashed lines. (a) Hyperbolic cosine-derived basis with increasing number of
basis functions indicated by darker shades of green (dark gray). (b) Piecewise
constant basis with increasing number of basis functions indicated by darker
shades of blue (light gray). The case $N=1$ corresponds to naive averaging of
the bulk dynamics and predicts no pattern for $h>1$, indicated by the striped
region.

Figure 8b shows the quantification of mean peak distances using a piecewise
constant basis for projecting the bulk dynamics [Fig. 7a,d,e]. The case $N=1$
corresponds to a naive averaging of the bulk dynamics and deviates
significantly from the reference value $\mu_{\Lambda}(h)$ already for small
bulk heights $h<\ell$. For $h\geq\ell$, the naive projection predicts no
patterns at all (shaded region) and is insufficient for capturing the
bulkd–boundary coupled reaction–diffusion dynamics. The predictions can again
be improved by introducing additional basis functions, where already $N=2$
basis functions approximately capture the characteristic shape of the
reference $\mu_{\Lambda}(h)$ and additional improvements are obtained for
$N=3$.

Importantly, even with just a single basis function meaningful approximations
of the actual dynamics can be achieved for an adequate choice of the basis.
This becomes conspicuous by comparing the mean peak distances
$\mu_{\Lambda}(h)$ in the case $N=1$ for a $\cosh$-derived basis and the
piecewise constant basis in Fig. 8 (light green and blue lines). With only one
basis function, the naive averaging cannot be used to approximate the system’s
dynamics, whereas for the hyperbolic cosine basis a single basis function is
sufficient to capture the geometry-dependent characteristics of the system.
This improvement is achieved by leveraging the a priori knowledge about the
system’s sensitivity to the bulk geometry and respecting this sensitivity in
the projection step.

Figure 9: (a) Schematic visualization of two distinct types of deformation:
top right, deformation by compression of the domain along the tangential
direction of the membrane, indicated by smaller distance between black
stripes. Bottom right, deformation by bending the membrane with curvature $K$
while keeping the metric on the membrane constant, indicated by equidistant
black stripes. Inset, definition of the bulk–boundary ratio
$\mathrm{d}V/\mathrm{d}A$. (b) Bulk–boundary ratio for varying bulk heights
and membrane curvatures. For studying the projection method in deforming
geometries, we consider a geometry that oscillates between flat (black dot)
and outward bent (open dot). Additional characteristic geometries are
indicated by gray dots. For large bulk heights and strong inward curvature,
the bulk–boundary ratio is ill-defined due to self-intersection of the bulk
domain (white area).

###  III.3 Curved and dynamic geometry

So far, we discussed the accuracy of the projection ansatz for a flat and
static geometry. How does this method perform in a dynamically deforming
geometry? In general, two types of deformations can be distinguished:
deformations that lead to a (local) in-plane compression or stretching of the
boundary domain while keeping the shape of the entire domain unchanged, and
deformations that keep the metric of the membrane invariant but cause out-of-
plane shape changes [Fig. 9a]. In the former case, the membrane is directly
affected by the deformation, making it difficult to disentangle the effects of
the bulk deformation from the membrane deformation. Instead, we here focus on
the latter case where the reaction–diffusion dynamics on the membrane are not
directly affected by the geometry deformation but only indirectly via the
coupling to the dynamic bulk. This allows to compare how well a certain basis
choice can capture bulk deformations.

To study such deforming geometries, we apply the projection method presented
in Section II.2 to the toy model in Eq. (12) that shows height-dependent
interrupted coarsening. For this, we choose a geometry parametrization
$\bm{r}(s,z;t)$ that leaves the metric on the boundary invariant (i.e.,
$\left|\partial_{s}\bm{r}(s,0;t)\right|=1$) but has a spatio-temporally
varying membrane curvature $K(s,t)$ [Fig. 9b, Appendix E]. Due to this dynamic
curvature, the local ratio of bulk volume to membrane area and thereby the
magnitude of the bulk field $c(\bm{r},t)$ now also vary in space and time
[57]. Note that in the static and flat case discussed above this ratio was
only affected by the bulk height $h$, since in this case the curvature is
$K=0$ at all times [Fig. 9b].

Figure 10: (a) Sample kymograph of the membrane concentration $m(s,t)$ for the
toy model system in Eq. (12) in a periodically deforming domain with membrane
curvature $K(s,t)=K_{\text{max}}/2\cdot(1{-}\cos(2\pi\,t/T))$ as indicated by
the sketches above the kymograph. The length of the membrane in all
deformation states is $L=10$, the bulk height is $h=1$. At strongly bent
states, the bulkd–boundary ratio is large and the ensuing low concentration of
the bulk field $c$ leads to fewer peaks compared to the flat state. (b) Mean
peak distance $\mu_{\Lambda}(t)$ for the reference system (black) and for a
projection using a polynomial basis $v_{k}=\\{1,(h{-}z)^{2}\\}$ (green/dark
gray), averaged over five samples with different initial conditions each with
system size $L=50$. For the reference system, results from individual samples
are indicated in gray. The deformation is identical to the one shown in panel
(a). (c) Mean peak distance $\mu_{\Lambda}(t)$ for the reference system
(black) and for piecewise constant basis with $N=2$ (light blue/light gray)
and $N=3$ (dark blue/gray). For $N=1$, no pattern forms for this bulk height
(not shown).

Similar to the static case, where the bulk height affects the mean peak
distance $\mu_{\Lambda}(h)$, the interrupted coarsening in a dynamic geometry
depends on the bulk–boundary ratio $\mathrm{d}V/\mathrm{d}A$ and thus on the
membrane curvature $K$ as shown in Fig. 9b. Modulating the membrane curvature
over time thus provides a straightforward method for probing the ability of
the projection method to capture effects of the dynamic geometry by comparing
the mean peak distance $\mu_{\Lambda}\bigl{(}K(t)\bigr{)}$. For this, we vary
the curvature periodically and homogeneously,

$K(s,t)=\frac{K_{\text{max}}}{2}\bigl{(}1-\cos(2\pi\,t/T)\bigr{)}$ (17)

with periodicity $T=10^{6}$ and maximum membrane curvature
$K_{\text{max}}=-1$. Figure 10a shows a typical kymograph of the membrane
concentration $m(s,t)$ during multiple deformation cycles. When the
deformation is strongest, the bulk–boundary ratio is maximal (equivalent to an
increased bulk height $h$, Fig. 8) and the mean distance between peaks is
large. As the geometry returns to the flat shape, the pattern adapts by
spawning additional peaks and thereby decreasing the mean peak distance
$\mu_{\Lambda}$.

In Fig. 10b,c we compare the periodic coarsening and spawning of peaks for
various basis choices, using a bulk height $h=1$ and maximum curvature
$K_{\text{max}}=-1$. Similar to the flat case, the naive projection using a
single basis function $v_{k}(z)=\\{1\\}$ fails to capture the height-dependent
dynamics and predicts no patterns at any deformation state at this bulk height
(not shown). To reproduce the reference dynamics, at least two basis functions
are required. Note that a hyperbolic cosine basis yields non-algebraic
contributions for the projected PDEs in deforming geometries and is therefore
not suitable for numerical implementation. Instead, we achieve qualitatively
matching results by expanding the hyperbolic cosine to lowest order, in this
case using a polynomial basis with $v_{k}=\\{1,(h{-}z)^{2}\\}$ [Fig. 10b]. For
sufficiently small deformations, the polynomial basis captures the reference
results with high accuracy, but considerable deviations are observed at strong
deformations. As an alternative, a piecewise constant basis with two or more
basis functions may be used [Fig. 10c]. Interestingly, the piecwise constant
basis captures the pattern statistics more robustly across the entire
deformation cycle than the polynomial basis even for only $N=2$ basis
functions. The reason for this lies in the polynomial basis being derived from
the steady state concentration profile in a _flat_ geometry. For strong
deformations, this steady state profile may change its shape significantly so
that the polynomial basis does not capture the relevant bulk gradients
anymore.

##  IV Discussion and Conclusion

In conclusion, we have presented a method for using a priori knowledge about
the physical characteristics of a bulk–boundary reaction–diffusion system to
reduce the complexity of the underlying model and thereby simplify the
numerical implementation of the system. Specifically, we proposed a systematic
approach to project the field dynamics in a volume onto the surface while
respecting the (possibly deforming) geometry of the volume and accounting for
spatial gradients normal to the surface. We showed that this method yields
more accurate approximations of the actual dynamics than a naive projection
based on averaging the field dynamics over the volume. At the same time, the
projection approach effectively lowers the dimension of the problem by one and
thus reduces the computational cost of numerically solving the
reaction–diffusion dynamics considerably (approximately 5-fold compared to
direct FEM implementations). As a proof of concept, we applied our method to a
generic model showing interrupted coarsening, where the statistical properties
of the coarsening dynamics are sensitive to the geometry of the volume. We
evaluated the accuracy of multiple projection methods in both static and
deforming geometries and found that for linear bulk reactions a set of two
basis functions already leads to good approximations of the actual dynamics
for all bulk heights. Using three or more basis functions provides minor
improvements that in general do not outweigh the increase in computation time
due to the additional fields.

We emphasize that the projection of the bulk dynamics onto the membrane as
presented in this article is only an approximation of the actual dynamics. In
contrast to finite element methods, where accuracy can often be increased
simply by refining the mesh [53], the accuracy of the projection strongly
depends on the basis choice. Leveraging the full potential of the projection
method therefore requires to make use of a priori knowledge about the physical
properties of the system at hand, such as the system’s steady state profiles
[54]. A physics-informed projection can, by virtue of the dimensionality
reduction, simplify the analytical and numerical assessment of pattern-forming
system. However, we stress that no significant benefits compared to FEM
implementations are to be expected for uninformed or arbitrary basis choices.

In our derivations of the projection method in deforming geometries we relied
on the geometry to be parametrizable and to have a constant height (as
measured from the membrane). As a natural extension to our results, it would
be interesting to release these constraints. In particular, the projection
method could be coupled to phase–field implementations of reaction–diffusion
dynamics on deforming surfaces, which have been extensively studied in the
past [58, 50, 46]. For this, a key challenge will be to implement a spatio-
temporally varying support of the basis functions to account for variable
volume height. Other extensions to our approach could include hydrodynamic
coupling of the field dynamics to the deforming geometry [22], which is an
important aspect of many intracellular reaction–diffusion systems [59, 60, 61,
62, 63], or basis functions that depend on time or are non-local in time, as
recently done for a phase-separating system with bulk–boundary coupling [64].

Beyond the modeling of reaction–diffusion dynamics in parametrized deforming
geometries, we expect our results to have valuable applications for systems
that reciprocally couple the geometry deformation to the field dynamics,
specifically for curvature-generating protein systems [65, 66]. Although the
projection method only yields an approximation of the true dynamics, we
believe that our approach can have significant advantages for large-scale
sampling of the model parameter space in simulations.

##  V Acknowledgements

We thank Tobias Roth and Henrik Weyer for stimulating discussions. T.B.
acknowledges support by the Joachim Herz Foundation. This work was funded by
the European Union (ERC, CellGeom, project number 101097810) and the Chan-
Zuckerberg Initiative (CZI).

##  Appendix A General derivation of projection methods

In this section, we provide a step-by-step derivation of the projection method
for arbitrarily deforming geometries. For additional information on finite
element methods, which are conceptually similar to the projection method
introduced here, we refer to pertinent textbooks [52, 53].

We start from a reaction–diffusion equation for a single field in a deforming
parametrizable geometry (specified by $\bm{r}(\bm{s},z,t)$, c.f. Fig. 4, with
the bulk extending in $z$-direction) with Robin boundary conditions at the
surface $\mathcal{S}$ (specified by
$\bm{r}_{\mathcal{S}}(\bm{s},t)=\bm{r}(\bm{s},0,t)$) and no-flux boundaries
elsewhere:

$\displaystyle\frac{1}{\sqrt{g}}\partial_{t}\bigl{(}\sqrt{g}\,c\bigr{)}$
$\displaystyle=D_{c}\,\Delta_{\text{LB}}c+f(c)\,,$ (18) $\displaystyle-
D_{c}\left.\bm{\hat{n}}(\bm{s},t)\cdot\bm{\nabla}c\right|_{\mathcal{S}}$
$\displaystyle=f_{\mathcal{S}}(\left.c\right|_{\mathcal{S}};\bm{s},t)\,,$ (19)
$\displaystyle\text{with}\qquad\Delta_{\text{LB}}c$
$\displaystyle=\frac{1}{\sqrt{g}}\partial_{i}\bigl{(}\sqrt{g}g^{ij}\partial_{j}c\bigr{)}\,.$
(20)

The functions $f$ and $f_{\mathcal{S}}$ denote the bulk reactions and the
bulkd–boundary coupling, respectively, and are both assumed to be polynomial
in $c$. We constrain ourselves to systems with constant bulk height $h$
measured in the normal direction at each point of the membrane. Consequently,
the $z$-direction is chosen to be normal to the membrane. This allows to write
the metric tensor in a block-diagonal form,

$g_{ij}=\begin{pmatrix}g_{\mathcal{S},i^{\prime}j^{\prime}}&0\\\
0&1\end{pmatrix}\,,$ (21)

where the indices $i,j$ run over all coordinates $\\{\bm{s},z\\}$ and the
indices $i^{\prime},j^{\prime}$ run only over the surface coordinates
$\\{\bm{s}\\}$ so that $g_{\mathcal{S},i^{\prime}j^{\prime}}$ is the metric
along the coordinates $\bm{s}$ of the reactive surface $\mathcal{S}$. In this
parametrization, the metric determinant $\det(g_{ij})$ is identical to the
determinant of the sub-metric,
$g\equiv\det(g_{\mathcal{S},i^{\prime}j^{\prime}})$.

Analogous to finite element methods [53], we start by transforming the PDE for
$c(\bm{s},z,t)$ to its weak form by introducing a test function $v(z)$ and an
inner product $\langle\cdot,\cdot\rangle^{ij}$ on the dynamic geometry along
the $z$-direction anchored at each point $\bm{s}$ on the reactive surface:

$\displaystyle\langle c(\bm{s},z,t),v(z)\rangle^{ij}$
$\displaystyle=\int_{0}^{z}\\!\\!\\!\mathrm{d}z\;\sqrt{g}\,g^{ij}\,c(\bm{s},z,t)\cdot
v(z)\,,$ (22) $\displaystyle\langle c(\bm{s},z,t),v(z)\rangle$
$\displaystyle=\int_{0}^{z}\\!\\!\\!\mathrm{d}z\;\sqrt{g}\,c(\bm{s},z,t)\cdot
v(z)\,,$ (23)

where we make use of the special form of the metric (21) to single out the
$z,z$ component of the inner product
$\langle\cdot,\cdot\rangle^{zz}=\langle\cdot,\cdot\rangle$ for easier
readability. With this, the individual terms of the partial differential
equation (18) in weak form can be expressed as

$\displaystyle\left\langle\frac{1}{\sqrt{g}}\partial_{t}\bigl{(}\sqrt{g}\,c\bigr{)},v\right\rangle$
$\displaystyle=\int_{0}^{z}\\!\\!\\!\mathrm{d}z\;\partial_{t}\bigl{(}\sqrt{g}\,c\bigr{)}\cdot
v=\left\langle\partial_{t}c,v\right\rangle+\left\langle\frac{1}{\sqrt{g}}c\cdot\partial_{t}\sqrt{g},v\right\rangle\,,$
(24) $\displaystyle\left\langle D_{c}\,\Delta_{\text{LB}}c,v\right\rangle$
$\displaystyle=D_{c}\int_{0}^{z}\\!\\!\\!\mathrm{d}z\;\left[\partial_{i^{\prime}}\left(\sqrt{g}\,g^{i^{\prime}j^{\prime}}\,\partial_{j^{\prime}}c\right)\cdot
v+\partial_{z}\left(\sqrt{g}\,\partial_{z}c\right)\cdot v\right]$
$\displaystyle=D_{c}\,\partial_{i^{\prime}}\left\langle\partial_{j^{\prime}}c,v\right\rangle^{i^{\prime}j^{\prime}}+D_{c}\left[\sqrt{g}\,\partial_{z}c\cdot
v\right]_{0}^{h}-D_{c}\langle\partial_{z}c,\partial_{z}v\rangle\,.$ (25)

Note that the index $i^{\prime}$ runs over all coordinates excluding the
$z$-coordinate, so that the derivative $\partial_{i^{\prime}}$ and the
integral over $z$ can be interchanged. In the next step, we apply a Galerkin
projection [53]: instead of solving the equations for arbitrary fields
$c(\bm{s},z,t)\in L^{2}(\mathcal{B})$, we only aim to find solutions for a
subspace $V\subset L^{2}(\mathcal{S})\times H^{1}([0,h])$. We choose this
subspace such that all functions $u(\bm{s},z,t)\in V$ can be expressed as

$u(\bm{s},z,t)=u_{k}(\bm{s},t)\,v_{k}(z)\,,$ (26)

where $v_{k}(z)$ is a set of $N$ basis functions in the Sobelov space
$H^{1}([0,h])$ with $k\in\\{1,\ldots,N{-}1\\}$, and we use Einstein summation
over double indices. To highlight the difference between indices denoting
spatial coordinates ($i,j$) and indices denoting elements of the basis
functions ($k,l,m$), we use co-/contravariant notation for the spatial
coordinates but consistently lower indices for the basis functions. For all
elements of the subspace $V$, the weak form of the differential equations can
be rewritten further via this separation ansatz (26):

$\displaystyle\left\langle\frac{1}{\sqrt{g}}\partial_{t}\bigl{(}\sqrt{g}\,c\bigr{)},v\right\rangle$
$\displaystyle=(\partial_{t}u_{k})\langle
v_{k},v\rangle+u_{k}\,\partial_{t}\langle v_{k},v\rangle\,,$ (27)
$\displaystyle\left\langle D_{c}\,\Delta_{\text{LB}}c,v\right\rangle$
$\displaystyle=D_{c}\,\partial_{i^{\prime}}\left[\bigl{(}\partial_{j^{\prime}}u_{k}\bigr{)}\,\langle
v_{k},v\rangle^{i^{\prime}j^{\prime}}\right]+$
$\displaystyle\qquad+D_{c}\left[\sqrt{g}\,\partial_{z}c\cdot
v\right]_{0}^{h}+$ $\displaystyle\qquad-
D_{c}\,u_{k}\langle\partial_{z}v_{k},\partial_{z}v\rangle\,.$ (28)

For the boundary term $[\sqrt{g}\,\partial_{z}c\cdot v]_{0}^{h}$, one may at
this point insert the boundary conditions $f_{\mathcal{S}}$ (for $h=0$) and
the no-flux boundary condition (for $h=z$). In the Galerkin projection, the
polynomial bulk reaction term $f(c)$ takes a special form. For the first and
second order, for example, the corresponding expressions are

$\displaystyle\langle c,v\rangle$ $\displaystyle=u_{k}\,\langle
v_{k},v\rangle\,,$ (29) $\displaystyle\langle c^{2},v\rangle$
$\displaystyle=u_{k}\,u_{l}\int_{0}^{z}\\!\\!\\!\mathrm{d}z\;\sqrt{g}\,v_{k}\,v_{l}\,v\,,$
(30)

and similarly for higher orders. Since the weak forms Eqs. (27)-(30) hold for
all test functions $v\in V$, one may conveniently choose the basis functions
$v_{k}$ to obtain PDEs for the coefficients $u_{k}(\bm{s},t)$. Introducing a
generalized Gram tensor $G_{kl}^{ij}$, a multilinear variant of the Gram
tensor $G^{ij}_{k_{1}\cdots k_{n}}$, and a generalized stiffness tensor
$A_{kl}^{ij}$ as

$\displaystyle G_{kl}^{ij}$
$\displaystyle=\int_{0}^{z}\\!\\!\\!\mathrm{d}z\;\sqrt{g}\,g^{ij}\,v_{k}\,v_{l}=\langle
v_{k},v_{l}\rangle^{ij}\,,$ (31) $\displaystyle G_{k_{1}\cdots k_{n}}^{ij}$
$\displaystyle=\int_{0}^{z}\\!\\!\\!\mathrm{d}z\;\sqrt{g}\,g^{ij}\,v_{k_{1}}\cdots
v_{k_{n}}\,,$ (32) $\displaystyle A_{kl}^{ij}$
$\displaystyle=\int_{0}^{z}\\!\\!\\!\mathrm{d}z\;\sqrt{g}\,g^{ij}\,(\partial_{z}v_{k})(\partial_{z}v_{l})$
$\displaystyle=\langle\partial_{z}v_{k},\partial_{z}v_{l}\rangle^{ij}\,,$ (33)

and abbreviating the entry $G^{zz}_{kl}=\langle v_{k},v_{l}\rangle\equiv
G_{kl}$ (similarly for $A_{kl}$ and $G_{k_{1}\cdots k_{n}}$) yields a concise
formulation for the dynamics of $u_{k}(\bm{s},t)$:

$(\partial_{t}u_{k})\,G_{kl}+\underbrace{u_{k}\,\partial_{t}G_{kl}\vphantom{\sum_{n}}}_{g_{\text{Geom}}}=\underbrace{D_{c}\,\partial_{i^{\prime}}\left[(\partial_{j^{\prime}}u_{k})G^{i^{\prime}j^{\prime}}_{kl}\right]\vphantom{\sum_{n}}}_{g_{\mathcal{S}\text{-Diff}}}-\underbrace{\left.\sqrt{g}\,f_{\mathcal{S}}(u_{k}\,v_{k})\cdot
v_{l}\right|_{0}\vphantom{\sum_{n}}}_{g_{\text{BC}}}-\underbrace{D_{c}\,u_{k}\,A_{kl}\vphantom{\sum_{n}}}_{g_{z\text{-Diff}}}+\underbrace{\sum_{n}\tfrac{1}{n!}\tfrac{\partial^{n}\\!f}{\partial
c^{n}}u_{k}\,u_{l_{2}}\cdots u_{l_{n}}\,G_{kl_{2}\cdots
l_{n}}}_{g_{\text{React}}}\,.$ (34)

In the case of purely linear bulk reactions, the reaction term reduces to
$g_{\text{React}}=\lambda\,u_{k}\,G_{kl}$ with linear reaction rate $\lambda$,
as applied in the toy model in Section III.

Importantly, for a parametrizable geometry, the generalized Gram tensor and
stiffness tensor can be calculated independently of the actual
reaction–diffusion dynamics. In particular, they can be calculated prior to
solving the partial differential equations for the coefficients $u_{k}$. The
task of solving the bulk dynamics is therefore largely shifted to calculating
the Gram and stiffness tensors, which greatly reduces the number of degrees of
freedom required to numerically solving the PDEs. It should be emphasized that
this is the primary distinguishing aspect compared to finite element methods:
in FEM implementations, the Galerkin projection is performed for each mesh
point in the bulk and on the boundary [53]. For a mesh size $a$ and system
width $L$ and height $h$ in $d$ dimensions, this results in a matrix of size
$\mathcal{O}\bigl{(}(L/a)^{d-1}\,(h/a)\bigr{)}$ that needs to be diagonalized
in each step. In contrast, by projecting the bulk dynamics onto the membrane
using $N$ basis functions the computational complexity is reduced to
diagonalizing a matrix of size $\mathcal{O}\bigl{(}N\,(L/a)^{d{-}1}\bigr{)}$,
where FEM (or other standard methods) can be used to solve the dynamics on the
membrane. Specifically, the projection method changes the size of the relevant
matrices by a factor $N\cdot a/h$, so that the computational benefit is
largest for low mesh sizes $a$ and a small number of basis functions $N$.

This improvement is possible by explicitly allowing highly nonlinear basis
functions in the projection approach and by adapting the basis functions
choice to the physical problem that is to be solved. For FEM, in contrast, it
is often convenient to choose basis functions $v_{k}$ with small support
centered around the mesh points, as this renders the Gram matrix $G$ as well
as other important matrices sparse and thereby speeds up numerical
applications. The projection method abolishes this need for basis functions
with small support in the bulk, since the Gram and stiffness matrices for the
bulk dynamics only need to be calculated once before solving the system
dynamics (however, Gram and stiffness matrix for the boundary dynamics still
need to diagonalized at each time step). When choosing the basis for the
projection method, one is therefore not constrained by the small support of
the bulk basis functions, and instead the basis can be tuned to the physical
properties of the system. A good choice of basis functions is therefore
crucial for leveraging the full potential of this dimensionality reduction
approach.

##  Appendix B Basis of orthogonal step functions

In this section, we discuss how the stiffness matrix $A_{kl}$ needs to be
adjusted when choosing a basis consisting of step functions. For this,
consider a one-dimensional bulk $[0,h]$ subdivided into $N$ segments at
$\\{0,z_{1},z_{2},\ldots,z_{N-1},h\\}$ so that each segment has size
$h_{k}=z_{k{+}1}-z_{k}$. The corresponding basis functions are

$v_{k}=\begin{cases}1&z_{k}\leq z\leq z_{k+1}\,,\\\ 0&\text{else}.\end{cases}$
(35)

This basis is orthogonal by design and thus the Gram matrix takes the form
$G=\text{diag}(\\{h_{0},\ldots,h_{N{-}1}\\})$. However, to construct the
stiffness matrix $A_{kl}$ it is necessary to find derivatives of the basis
functions or at least derivatives in the weak sense [53], i.e., functions
$w_{k}(z)$ for which

$\int_{0}^{h}\\!\\!\mathrm{d}z\;v_{k}(z)\cdot(\partial_{z}\phi(z))=-\int_{0}^{h}\\!\\!\mathrm{d}z\;w_{k}(z)\cdot\phi(z)$
(36)

for all differentiable test functions $\phi(z)$. Since step functions are not
weakly differentiable [67], the stiffness matrix $A_{kl}$ as defined in Eq.
(33) is not applicable for a step function basis.

At first glance, one may try to circumvent this issue by approximating the
step functions using smooth analogues, for example a logistic sigmoid
$\sigma_{w}(z)=(1+e^{-z/w})^{-1}$, where the parameter $w$ quantifies the
width of the step. The basis elements may then be written as
$v_{k}=\sigma_{w}(z-z_{k})\,\sigma_{w}(z_{k{+}1}-z)$. This basis is not
orthogonal anymore, however the Gram basis remains to lowest order
$G=\text{diag}(\\{h_{0},\ldots,h_{N{-}1}\\})+\mathcal{O}(w)$. The stiffness
matrix, on the other hand, takes a tridiagonal form to lowest order,

$A=\frac{1}{6w}\begin{pmatrix}-2&1&0&\\\ 1&-2&1&\cdots\\\ 0&1&-2&\\\
&\vdots&&\ddots\end{pmatrix}+\mathcal{O}(w)\,.$ (37)

As a result, the diffusion term $g_{\text{Diff}}$ with such a stiffness matrix
enters the differential equation for the coefficients (e.g., Eq. (5)) with a
prefactor ${\sim}w^{-1}$, leading to diverging terms in the limit $w\to 0$.
This is because the true gradients $\partial_{z}c(z)$ are poorly approximated
by the step functions, where gradients are
$\partial_{z}u\approx(u_{k{+}1}-u_{k})/w$ in the $w$-neighborhood of the nodes
$z_{k}$ and zero elsewhere. To be able to use step functions (or their smooth
analogues) as a basis, the stiffness matrix needs to be reconstructed
artificially.

For this, it is necessary to define approximations of the gradient terms
$\partial_{z}\tilde{v}_{k}$ as stand-ins for the weak derivatives of the
actual basis $v_{k}$. As the simplest approach, one may emulate standard FEM
strategies and use linear interpolation between the separate layers [Fig. XX].
This yields piecewise constant functions for the approximate gradients:

$\partial_{z}\tilde{v}_{k}=\begin{cases}\frac{2}{h_{k{-}1}+h_{k}}&z_{k}-\frac{h_{k{-}1}}{2}\leq
z\leq z_{k}+\frac{h_{k}}{2}\\\
\frac{2}{h_{k}+h_{k{+}1}}&z_{k{+}1}-\frac{h_{k}}{2}\leq z\leq
z_{k{+}1}-\frac{h_{k+1}}{2}\\\ 0&\text{else}\,.\end{cases}$ (38)

The resulting reconstructed stiffness matrix $\bar{A}_{kl}$ then evaluates to
the tridiagonal form

$\displaystyle\bar{A}_{kk}$
$\displaystyle=\frac{2}{h_{k}+h_{k{+}1}}+\frac{2}{h_{k}+h_{k{-}1}}\,,$
$\displaystyle\bar{A}_{k,k{\pm}1}$
$\displaystyle=-\frac{2}{h_{k}+h_{k{\pm}1}}\,,$ (39)

and zero for all other entries. A priori knowledge about the system, e.g.,
about the steady state distribution, can be used to obtain better estimates
for the gradients. For all results presented in this paper we use the FEM-
inspired stiffness matrix as defined in Eq. (39) when working with the step
function basis.

##  Appendix C Numerical implementation

Parameter | Value |  Description
---|---|---
$D_{c}$ | $1$ |  Bulk diffusion
$D_{m}$ | $10^{-4}$ |  Membrane diffusion
$\lambda$ | $1$ |  Bulk degradation rate
$\epsilon$ | $10^{-2}$ |  Membrane degradation rate
$n_{0}$ | $1$ |  Average initial mass on membrane
$L$ | $100$ |  System width
$\ell=\sqrt{D_{c}/\lambda}$ | $1$ |  Characteristic length scale
$p$ | $1.1$ |  Effective membrane activation rate
Table 1: List of parameters used in simulations.

To test the projection method with the toy model as presented in Section III,
we performed numerical simulations using COMSOL Multiphysics v6.0 for the
different projection approaches (basis choices) and, for reference, for a non-
projected system, with simulation parameters as stated in Table 1. For the
reference simulations, the bulk dynamics were solved in a two-dimensional
domain coupled to a one-dimensional boundary accounting for the membrane
dynamics via reactive boundary conditions. The boundary domain was meshed
regularly with a mesh size of $0.05$. The bulk domain was subsequently meshed
using a triangular Delauney tesselation (irregular) combined with local
adaptive mesh refinement, where the mesh resolution was reduced at distances
far away from the membrane controlled by setting the software’s configuration
parameter Maximum element growth rate to $1.5$. Note that the adaptive mesh
refinement already assumes shallow gradients far away from the membrane. For
comparison, we also produced solutions for a bulk domain that was meshed in a
regular rectangular grid with mesh size $0.05$. For the projected systems, all
dynamics (membrane and projected bulk) were solved on a one-dimensional domain
corresponding to the boundary in the reference simulations. This domain was
meshed regularly with a mesh size of $0.05$.

Figure 11: Benchmark of solution times for the toy model in a flat geometry
for a system with explicit bulk dynamics (reference, black, with regular
rectangular mesh (squares) and local adaptive triangular mesh (triangles)),
for a projected system with basis functions derived from the hyperbolic cosine
steady state concentration profile (green/dark gray, diamonds), and for a
projected system with piecewise constant basis functions (blue/light gray,
circles). Horizontal displacement is added to ensure that data points are
distinguishable but does not correspond to different bulk heights $h$. Data
points represent the average value over five samples each. Relative standard
deviations ($\sigma_{t}/\mu_{t}$) are $<0.45$ (reference, $h=10$) and $<0.22$
for all other values. Lines are added to guide the eye.

In Fig. 11 we compare the solution times for different simulation methods
(reference and projections) for three different values of the system height
$h$. Each data point in Fig. 11 is the mean of five samples obtained from
different initial conditions and otherwise identical parameters. In the
reference simulations (both regular mesh and adaptive refined mesh, black),
solution times are consistently larger than ten seconds, and increase to
approximately one minute (adaptive refined mesh) and ten minutes (regular
mesh) at large system heights $h=10$ where the number of mesh points becomes
large. In contrast, for simulations using bulk projection with up to $N=3$
basis functions the solution time remained below 20 seconds, and reached
approximately five seconds for low system heights ($h=0.1$, where coarsening
is interrupted earlier and a steady state is reached faster, c.f. Fig. 6).
Importantly, the number of mesh points in the projected system is independent
of the bulk height in contrast to the reference system, and therefore the
computational speedup is most prominent for large bulk height and a regular
mesh. Comparing the reference simulations with adaptive refined mesh to the
simulations using $N=1$ basis functions derived from a hyperbolic cosine, a
speedup factor of ${>}3$ is achieved for $h=0.1$, and a factor ${>}4$ for
$h=1$ and $h=10$. For benchmarking the solution times, the simulations were
performed on a 64 core AMD Ryzen Threadripper processor at 2.70 GHz.

##  Appendix D Projection of mass-conserving systems

In reaction–diffusion systems, the diffusion part always conserves the total
protein mass $n=\int\\!\mathrm{d}V\,c(\bm{x},t)$, i.e., for a system with no-
flux boundaries at all boundaries and without any reaction terms defined by

$\partial_{t}c(\bm{x},t)=D_{c}\,\nabla^{2}c(\bm{x},t)$ (40)

the change in total mass is [3]

$\partial_{t}n=\int_{\mathcal{B}}\\!\\!\mathrm{d}V\;\partial_{t}c(\bm{x},t)=\left.\nabla
c(\bm{x},t)\right|_{\partial\mathcal{B}}=0\,.$ (41)

This property of mass-conserving diffusion should be preserved under the
projection method, in particular for systems where all reactions are also
mass-conserving. Under what circumstances does the projection method respect
mass conservation in the diffusion?

For this, it is sufficient to consider diffusion in a one-dimensional bulk
(i.e., no diffusion parallel to the membrane) with no-flux boundary conditions
at both boundaries. For such a system, the governing equations for the fields
$u_{k}(t)$ derived in Eq. (34) reduce to

$G_{kl}\,(\partial_{t}u_{k})=-D_{c}\,u_{k}\,A_{kl}\,.$ (42)

Without loss of generality, assume furthermore that the basis $\\{v_{k}\\}$ is
orthonormal so that the Gram matrix is a unit matrix $G_{kl}=\delta_{kl}$.
Since each component $k$ contributes to the total concentration by
$u_{k}\,v_{k}(z)$, the change in the total mass is then

$\partial_{t}n=\partial_{t}u_{k}\,\int\\!\mathrm{d}z\;v_{k}(z)=-D_{c}\left[\int\\!\mathrm{d}z\;v_{k}(z)\right]\cdot
A_{kl}\,u_{l}\,.$ (43)

The total mass is thus conserved if and only if

$A_{kl}\,\int\\!\mathrm{d}z\;v_{k}(z)=0\quad\forall\,l\,.$ (44)

This poses an additional constraint on the basis choice, complementing the
requirements to approximate the true distribution $c(z,t)$ and its gradient
$\partial_{z}c(z,t)$ well. We now show that this requirement is always
fulfilled if one of the basis functions is constant in $z$-direction. For
this, assume that $v_{0}(z)=1/\sqrt{h}$ in a bulk domain of height $h$. This
allows to rewrite the integral over a single basis function in terms of the
Gram matrix,

$\int\\!\mathrm{d}z\;v_{k}(z)=\sqrt{h}\int\\!\mathrm{d}z\;v_{k}(z)\,v_{0}(z)=\sqrt{h}\,G_{k0}\,.$
(45)

Inserting this into the requirement for mass conservation in Eq. (44) and
using the orthonormality of the basis $G_{kl}=\delta_{kl}$ yields

$\sqrt{h}A_{kl}\,G_{k0}=\sqrt{h}A_{0l}=0\,,$ (46)

which is always fulfilled for a basis that contains a constant basis function
(in this case $v_{0}$). This follows from the fact that
$A_{kl}=\langle\partial_{z}v_{k},\,\partial_{z}v_{l}\rangle$, and by choice of
the basis function $\partial_{z}v_{0}=0$.

Since the orthonormalization step only affects the representation of the
system but not the dynamics, the same line of argument holds for a basis where
only a linear combination of basis functions is constant, i.e., where

$\alpha_{k}\,\partial_{z}v_{k}(z)=0$ (47)

for real constants $\alpha_{k}$. In particular, this includes a set of
piecewise constant basis functions as presented in Sec. III.

##  Appendix E Parametrization of evolving geometries

In this section, we specify the parametrization used to evaluate the
projection method in a deforming geometry as discussed in Sec. III.3. Since we
use that the basis functions $v_{k}(z)$ are identical for all points
$\bm{r}(\bm{s},0)$ on the membrane $\mathcal{S}$, a necessary constraint for
the parametrizations is that they have constant height $h$ everywhere. Note
that here the height is defined as the distance from the membrane
$\mathcal{S}$ to the opposing no-flux boundary in the direction normal to
$\mathcal{S}$. This also implies that these two boundaries are always
parallel. In addition, deformations need to avoid intersections of the bulk
with itself, since in such cases the metric diverges. Therefore, the maximum
height of the bulk in a deforming geometry is constrained by the maximum
(principal) curvature, $h<1/K_{\text{max}}$. Furthermore, we here choose a
parametrization for which the metric on the membrane is constant and unity,
$\left.g_{ij}\right|_{z=0}=\delta_{ij}$, since this ensures that no stretching
or compression of the membrane can affect the coarsening dynamics and all
observed effects are due to the bulk dynamics.

For the examples in Sec. III.3, we examine a geometry consisting of an annular
segment with spatially constant but temporally varying curvature $K(t)$. This
geometry is parametrized by

$\bm{r}(s,z,t)=\frac{1}{K(t)}\begin{pmatrix}\bigl{(}1+z\,K(t)\bigr{)}\sin(s\,K(t))\\\
\bigl{(}1+z\,K(t)\bigr{)}\cos(s\,K(t))-1\end{pmatrix}\,,$ (48)

where we choose
$K(t)=K_{\text{max}}\cdot\frac{1}{2}\left((1-\cos(\frac{2\pi\,t}{T})\right)$
with $K_{\text{max}}<0$. This geometry is flat at $t=n\cdot 2\pi\,T$ for
$n\in\mathbb{N}$ and morphs into an arc with curvature $|K_{\text{max}}|$ at
$t=(2n+1)\cdot\pi\,T$ with a periodicity $T$. This parametrization ensures
that the surface is never inward bent and therefore no problematic self-
intersections can occur.

## References

  * Wolpert [1981] L. Wolpert, Philosophical Transactions of the Royal Society of London. B, Biological Sciences 295, 441 (1981).
  * Desai and Kapral [2009] R. C. Desai and R. Kapral,  _Dynamics of Self-Organized and Self-Assembled Structures_ (Cambridge University Press, 2009).
  * Burkart  _et al._ [2022] T. Burkart, M. C. Wigbers, L. Würthner, and E. Frey, Nature Reviews Physics 4, 511 (2022).
  * Raskin and Boer [1999] D. M. Raskin and P. A. J. d. Boer, Proceedings of the National Academy of Sciences 96, 4971 (1999).
  * Ramm  _et al._ [2019] B. Ramm, T. Heermann, and P. Schwille, Cellular and Molecular Life Sciences 76, 4245 (2019).
  * Chiou  _et al._ [2016] J.-G. Chiou, M. K. Balasubramanian, and D. J. Lew, Annual Review of Cell and Developmental Biology 33, 1 (2016).
  * Motegi and Seydoux [2013] F. Motegi and G. Seydoux, Philosophical Transactions of the Royal Society B: Biological Sciences 368, 20130010 (2013).
  * Walani  _et al._ [2014] N. Walani, J. Torres, and A. Agrawal, Physical Review E 89, 062715 (2014).
  * Rogez  _et al._ [2019] B. Rogez, L. Würthner, A. B. Petrova, F. B. Zierhut, D. Saczko-Brack, M.-A. Huergo, C. Batters, E. Frey, and C. Veigel, Nature Communications 10, 3305 (2019).
  * Mahapatra  _et al._ [2021] A. Mahapatra, D. Saintillan, and P. Rangamani, Soft Matter 17, 8373 (2021).
  * Roux  _et al._ [2021] A.-L. L. Roux, C. Tozzi, N. Walani, X. Quiroga, D. Zalvidea, X. Trepat, M. Staykova, M. Arroyo, and P. Roca-Cusachs, Nature Communications 12, 6550 (2021).
  * Würthner  _et al._ [2023] L. Würthner, A. Goychuk, and E. Frey, Physical Review E 108, 014404 (2023), 2205.02820 .
  * Hubatsch  _et al._ [2019] L. Hubatsch, F. Peglion, J. D. Reich, N. T. L. Rodrigues, N. Hirani, R. Illukkumbura, and N. W. Goehring, Nature Physics 15, 1078 (2019).
  * Marée  _et al._ [2012] A. F. M. Marée, V. A. Grieneisen, and L. Edelstein-Keshet, PLoS Computational Biology 8, e1002402 (2012).
  * Camley  _et al._ [2017] B. A. Camley, Y. Zhao, B. Li, H. Levine, and W.-J. Rappel, Physical Review E 95, 012401 (2017), 1609.01764 .
  * Fu  _et al._ [2023] M. Fu, T. Burkart, I. Maryshev, H. G. Franquelim, A. Merino-Salomón, M. Reverte-López, E. Frey, and P. Schwille, Nature Physics 19, 1211 (2023).
  * Miller  _et al._ [2018] P. W. Miller, N. Stoop, and J. Dunkel, Physical Review Letters 120, 268001 (2018), 1710.02247 .
  * Litschel  _et al._ [2018] T. Litschel, B. Ramm, R. Maas, M. Heymann, and P. Schwille, Angewandte Chemie International Edition 57, 16286 (2018).
  * Tamemoto and Noguchi [2020] N. Tamemoto and H. Noguchi, Scientific Reports 10, 19582 (2020), 2007.10606 .
  * Hörning and Shibata [2019] M. Hörning and T. Shibata, Biophysical Journal 116, 372 (2019).
  * Brinkmann  _et al._ [2018] F. Brinkmann, M. Mercker, T. Richter, and A. Marciniak-Czochra, PLOS Computational Biology 14, e1006259 (2018).
  * Mietke  _et al._ [2018] A. Mietke, F. Jülicher, and I. F. Sbalzarini, Proceedings of the National Academy of Sciences 116, 29 (2018).
  * Meacci  _et al._ [2006] G. Meacci, J. Ries, E. Fischer-Friedrich, N. Kahya, P. Schwille, and K. Kruse, Physical Biology 3, 255 (2006), q-bio/0701046 .
  * Dziuk and Elliott [2007] G. Dziuk and C. M. Elliott, IMA Journal of Numerical Analysis 27, 262 (2007).
  * Madzvamuse  _et al._ [2015] A. Madzvamuse, A. H. W. Chung, and C. Venkataraman, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471, 20140546 (2015), 1501.05769 .
  * Chen and Ling [2020] M. Chen and L. Ling, Journal of Computational Physics 405, 109166 (2020), 1906.01370 .
  * Bergdorf  _et al._ [2010] M. Bergdorf, I. F. Sbalzarini, and P. Koumoutsakos, Journal of Mathematical Biology 61, 649 (2010).
  * Hansbo  _et al._ [2016] P. Hansbo, M. G. Larson, and S. Zahedi, Computer Methods in Applied Mechanics and Engineering 307, 96 (2016), 1502.07142 .
  * Suchde and Kuhnert [2019] P. Suchde and J. Kuhnert, Journal of Computational Physics 395, 38 (2019), 1902.08107 .
  * Wendland and Künemund [2020] H. Wendland and J. Künemund, Advances in Computational Mathematics 46, 64 (2020).
  * Lervåg and Lowengrub [2015] K. Y. Lervåg and J. Lowengrub, Communications in Mathematical Sciences 13, 1473 (2015), 1407.7480 .
  * Loose  _et al._ [2008] M. Loose, E. Fischer-Friedrich, J. Ries, K. Kruse, and P. Schwille, Science 320, 789 (2008).
  * Bement  _et al._ [2015] W. M. Bement, M. Leda, A. M. Moe, A. M. Kita, M. E. Larson, A. E. Golding, C. Pfeuti, K.-C. Su, A. L. Miller, A. B. Goryachev, and G. v. Dassow, Nature Cell Biology 17, 1471 (2015).
  * Paquin-Lefebvre  _et al._ [2020] F. Paquin-Lefebvre, B. Xu, K. L. DiPietro, A. E. Lindsay, and A. Jilkine, Journal of Theoretical Biology 497, 110242 (2020), 1910.10203 .
  * Halatek  _et al._ [2018] J. Halatek, F. Brauns, and E. Frey, Philosophical Transactions of the Royal Society B: Biological Sciences 373, 20170107 (2018), 1802.07169 .
  * Brauns  _et al._ [2021a] F. Brauns, G. Pawlik, J. Halatek, J. Kerssemakers, E. Frey, and C. Dekker, Nature Communications 12, 3312 (2021a).
  * Würthner  _et al._ [2022] L. Würthner, F. Brauns, G. Pawlik, J. Halatek, J. Kerssemakers, C. Dekker, and E. Frey, Proceedings of the National Academy of Sciences 119, e2206888119 (2022), 2111.12043 .
  * Levine and Rappel [2005] H. Levine and W.-J. Rappel, Physical Review E 72, 061912 (2005).
  * Li  _et al._ [2009] X. Li, J. Lowengrub, A. Rätz, and A. Voigt, Communications in Mathematical Sciences 7, 81 (2009).
  * Teigen  _et al._ [2009] K. E. Teigen, X. Li, J. Lowengrub, F. Wang, and A. Voigt, Communications in Mathematical Sciences 7, 1009 (2009).
  * Yu  _et al._ [2012] H.-C. Yu, H.-Y. Chen, and K. Thornton, Modelling and Simulation in Materials Science and Engineering 20, 075008 (2012), 1107.5341 .
  * Marth and Voigt [2014] W. Marth and A. Voigt, Journal of Mathematical Biology 69, 91 (2014).
  * Rätz and Röger [2014] A. Rätz and M. Röger, Nonlinearity 27, 1805 (2014), 1305.6172 .
  * Abels  _et al._ [2015] H. Abels, K. F. Lam, and B. Stinner, SIAM Journal on Mathematical Analysis 47, 3687 (2015).
  * Moure and Gomez [2020] A. Moure and H. Gomez, Soft Matter 16, 5177 (2020).
  * Valizadeh and Rabczuk [2019] N. Valizadeh and T. Rabczuk, Computer Methods in Applied Mechanics and Engineering 351, 599 (2019).
  * Biben  _et al._ [2005] T. Biben, K. Kassner, and C. Misbah, Physical Review E 72, 041921 (2005).
  * Rätz and Voigt [2006] A. Rätz and A. Voigt, Communications in Mathematical Sciences 4, 575 (2006).
  * Elliott  _et al._ [2011] C. M. Elliott, B. Stinner, V. Styles, and R. Welford, IMA Journal of Numerical Analysis 31, 786 (2011).
  * Poulsen and Voorhees [2018] S. O. Poulsen and P. W. Voorhees, International Journal of Computational Methods 15, 1850014 (2018).
  * Yu  _et al._ [2020] F. Yu, Z. Guo, and J. Lowengrub, Journal of Computational Physics 406, 109174 (2020).
  * Zienkiewicz  _et al._ [2005] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhou,  _The finite element method: its basis and fundamentals_ (Elsevier, 2005).
  * Brenner and Scott [2008] S. C. Brenner and L. R. Scott, Texts in Applied Mathematics 10.1007/978-0-387-75934-0 (2008).
  * Frey and Brauns [2022] E. Frey and F. Brauns, in  _Active Matter and Nonequilibrium Statistical Physics_, Lecture Notes of the Les Houches Summer School, Vol. 112, edited by J. Tailleur (Oxford Academic, Oxford, 2022) pp. 347–445.
  * Brauns  _et al._ [2020] F. Brauns, J. Halatek, and E. Frey, Physical Review X 10, 041036 (2020), 1812.08684 .
  * Brauns  _et al._ [2021b] F. Brauns, H. Weyer, J. Halatek, J. Yoon, and E. Frey, Physical Review Letters 126, 104101 (2021b), 2005.01495 .
  * Thalmeier  _et al._ [2016] D. Thalmeier, J. Halatek, and E. Frey, Proceedings of the National Academy of Sciences 113, 548 (2016).
  * Ayton  _et al._ [2005] G. S. Ayton, J. L. McWhirter, P. McMurtry, and G. A. Voth, Biophysical Journal 88, 3855 (2005).
  * Goehring  _et al._ [2011] N. W. Goehring, P. K. Trong, J. S. Bois, D. Chowdhury, E. M. Nicola, A. A. Hyman, and S. W. Grill, Science 334, 1137 (2011).
  * Bois  _et al._ [2011] J. S. Bois, F. Jülicher, and S. W. Grill, Physical Review Letters 106, 028103 (2011).
  * Nishikawa  _et al._ [2017] M. Nishikawa, S. R. Naganathan, F. Jülicher, and S. W. Grill, eLife 6, e19595 (2017).
  * Klughammer  _et al._ [2018] N. Klughammer, J. Bischof, N. D. Schnellbächer, A. Callegari, P. Lénárt, and U. S. Schwarz, PLoS computational biology 14, e1006588 (2018).
  * Illukkumbura  _et al._ [2020] R. Illukkumbura, T. Bland, and N. W. Goehring, Current Opinion in Cell Biology 62, 123 (2020).
  * Caballero  _et al._ [2023] N. Caballero, K. Kruse, and T. Giamarchi, Physical Review E 108, L012801 (2023), 2205.03306 .
  * Alimohamadi and Rangamani [2018] H. Alimohamadi and P. Rangamani, Biomolecules 8, 120 (2018).
  * Tozzi  _et al._ [2019] C. Tozzi, N. Walani, and M. Arroyo, New Journal of Physics 21, 093004 (2019).
  * Adams and Fournier [2003] R. A. Adams and J. J. F. Fournier,  _Sobolev Spaces_ , 2nd ed. (Elsevier, 2003).