{'fraction_non_alphanumeric': 0.03379961827450412, 'fraction_numerical': 0.011089367253750815, 'mean_word_length': 4.661605799480236, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Drift-Diffusion Dynamics and Phase Separation in Curved Cell Membranes and Dendritic Spines: Hybrid Discrete-Continuum Methods

We develop methods for investigating protein drift-diffusion dynamics in heterogeneous cell membranes and the roles played by geometry, diffusion, chemical kinetics, and phase separation. Our hybrid stochastic numerical methods combine discrete particle descriptions with continuum-level models for tracking the individual protein drift-diffusion dynamics when coupled to continuum fields. We show how our approaches can be used to investigate phenomena motivated by protein kinetics within dendritic spines. The spine geometry is hypothesized to play an important biological role regulating synaptic strength, protein kinetics, and self-assembly of clusters. We perform simulation studies for model spine geometries varying the neck size to investigate how phase-separation and protein organization is influenced by different shapes. We also show how our methods can be used to study the roles of geometry in reaction-diffusion systems including Turing instabilities. Our methods provide general approaches for investigating protein kinetics and drift-diffusion dynamics within curved membrane structures.


I. INTRODUCTION
In cellular biology, the morphological shapes of cell membranes play important roles in protein transport and kinetics. Cell membranes often take on shapes having characteristic geometries or topologies associated with biological function [1][2][3][4][5][6][7]. Membranes arising in cell biology consist of heterogeneous mixtures of lipids, proteins, and other small molecules [1]. The individual and collective dynamics of membrane associated molecules carry out diverse functions in cellular processes ranging from signaling to motility [1,[8][9][10][11][12]. Membranes are effectively two dimensional fluid-elastic structures resulting in processes that can be significantly different than their counter-parts occuring in bulk three dimensional fluids [13][14][15][16]. Investigating such cellular processes using computational simulation requires the ability to capture these effects and the geometric and topological contributions of curved membrane structures to protein driftdiffusion dynamics and kinetics.
We introduce computational methods based on a hybrid approach coupling discrete and continuum descriptions. For low concentration species, we track individual proteins as discrete particles. For other species, we track contributions using continuum fields. We circumvent many of the challenges of differential geometry and directly approximating surface PDEs by developing discrete localized models that capture geometric effects. In this way, behaviors of our model emerge on larger lengthscales in a manner capturing the relevant underlying physical phenomena, while avoiding some of the more common challenges associated with direct application of PDE discretizations and differential geometry. We also provide stochastic local models to account for discrete effects and other fluctuations.
Our methods provide new ways for handling hybrid discrete-continuum stochastic descriptions for protein drift-diffusion dynamics and kinetics within heterogeneous curved membranes. We focus particularly on the roles played by phase separation, geometry, and fluctuations arising from discrete number effects. We develop stochastic methods for bidirectional coupling between protein dynamics and evolving phase fields. For capturing geometric contributions for general shapes, we also introduce numerical approaches using the induced metric from the embedding space which allows for avoiding the need for potentially cumbersome calculations using explicit expressions from differential geometry. We also introduce methods for tracking proteins both at the level of continuum mean-field concentration fields and beyond mean-field theory at the level of discrete individual proteins with stochastic trajectories.
Our work is motivated by understanding protein interactions within cell curved membranes, such as neuronal dendritic spines. Dendritic spines are small ∼ 500nm structures that are attached along the larger shaft of the dendrites of neurons [45,46]. These are critical structures in the brain mediating the input of synaptic communications between neurons. Dendritic spines have unique morphologies that physiologically change in molecular composition, size, and shape to modulate the synaptic strength between neurons as part of long-term learning and memory [46][47][48][49]. The important connection between membrane geometry and synaptic strength is an active area of current theoretical and experimental research [50][51][52][53][54][55][56][57][58][59][60][61][62]. Recent advances in microscopy and single-particle tracking techniques are providing some indicators of the underlying processes regulating protein transport and kinetics within spines [47,[63][64][65][66][67][68][69].
A hypothesis which we explore with our methods is that protein complexes can locally nucleate liquid-liquid phase separations that couple to proteins to influence both diffusive transport and kinetics. Our work is motivated by the recent observations concerning the roles of SynGAP binding to PSD-95 resulting in complexes that phase separate playing a role in organizing proteins in dendritic spines in forming the post-synaptic density [70,71]. In our initial work presented here, we do not commit in models yet to specific proteins, but focus more generally on mechanisms by which protein diffusion and phase separation can drive cluster formation and augment reaction kinetics consistent with such observations. For exploring such hypotheses, we focus on the development of methods for quantitative biophysical models and computational simulations capable of investigating such effects. For this purpose, we develop spatial-temporal biophysical models using our introduced hybrid discretecontinuum approach to investigate the roles of geometry and phase separation in protein diffusion and kinetics.
We discuss details for our modeling approach for protein drift-diffusion dynamics in Section II and our simulation studies in Section III. We discuss the drift-diffusion dynamics of proteins coupled to phase fields and discretizations for curved membranes having general shapes in Sections II A, II B, II C, and II D. We investigate the accuracy of our numerical approximations in Section II E. We demonstrate how our methods can be used to compute first passage times and other statistics for membranes of general shape, in some cases without the need for costly Monte-Carlo sampling, in Section II F.
We perform simulations using our methods of discrete and continuum systems in III. We develop continuum field methods for reaction-diffusion systems with pattern forming Turing instabilities capturing the contributions of geometry in Section III A. We develop mechanistic models for dendritic spine shapes and perform simulation to investigate the roles of geometric effects, heterogeneities arising from phase separation, and discrete number effects in Section III B and Section III C. The introduced methods provide general approaches for investigating protein transport and kinetics within heterogeneous cell membranes. The surface is discretized using a triangulation with random walks between centriods xi and xj having jump rates Mij determined by the local geometry (top-left) which can be expressed as a Markov-Chain (bottom). The Markov-Chain discretization yields a backward equation for Ω allowing readily for computation of firstpassage times τ (FPTs) for reaching the boundary ∂Ω and other statistics without the need in some cases for Monte-Carlo sampling (top-right).

II. HYBRID DISCRETE PARTICLE-CONTINUUM APPROACH FOR CURVED SURFACES
We develop models using a hybrid discrete particle and continuum field approach for investigating transport and reactions within curved surfaces, see Figure 1. For low concentration chemical species, we introduce approaches for modeling at the level of tracking individual particles. For larger concentrations, we develop continuum field descriptions and methods for curved surfaces. We also account for the bi-directional coupling between the discrete particles and continuum fields.

A. Particle Drift-Diffusion Dynamics
We discuss a few results related to the drift-diffusion dynamics of particles useful in developing our methods. We then approximate the particle dynamics using a Markov-Chain process [72] with jump rates based on estimating the local geometry, see Figure 1.
The drift-diffusion dynamics of a particle immersed within a viscous fluid is given by the Langevin equation [73,74] where dX t = V t dt. This is to be interpreted as an Itô process [73,75]. The γ is the drag, U (x) is a potential energy, k B is Boltzmann's constant, and T is the temperature [74]. The dW t are increments of the Wiener process [73,75]. The diffusion coefficient is given by We remark that in some cases the protein diffusion can also be influenced by local protein-induced deformations of the membrane [16], crowding effects (sub-super diffusive and other regimes) [76,77], or cytoskeletal interactions [24,78]. Our methods could be combined with such models for local deformation and crowding by augmenting the Markov-Chain jump rates. In the present work, we treat the membrane as having on average an effectively constant shape. We focus on models in the normal diffusive regime with explicit modeling of protein species, binding partners, or obstacles. We treat additional contributions through effective diffusion coefficients or in the formulated protein interaction forces.
When m/γ 2 /D, where is the radius of the particle, the inertial contributions are negligible. In this regime, the Langevin equation can be reduced to the over-damped Smoluchowski equation This can be expressed in terms of probability densities ρ(x, t) for when X t = x. This satisfies the Fokker-Planck (FP) equation When U = 0, this has the well-known Green's Function for Euclidean space From the FP equation, the K(x , x; t) has the interpretation of the probability of a particle starting with X 0 = x and diffusing to location X t = x over the time duration t. We will use a related approach to determine our Markov-Chain jump rates. The FP equation has a steady-state ρ * with detailed balance if For a smooth U (x) with a sufficient growth rate as |x| → ∞, the FP equation has steady-state with detailedbalance for the distribution This is the Gibbs-Boltzmann distribution, and Z is the partition function normalizing this to be a probability density [74].

B. Markov-Chain Discretization for Particle Drift-Diffusion Dynamics on Curved Surfaces
We model the drift-diffusion dynamics of individual particles on curved surfaces using discrete Markov-Chains with the jump rates based on the local geometry. The surface is discretized into a triangulated mesh and each particle is tracked by the triangle which it occupies. We use that the surface metric is induced by the surrounding embedding space [79,80]. We use for the local jump rates To approximate the diffusion over the time-scale ∆t on the surface M, we use = √ 4D∆t. The x i ∈ M are the centers of the surface triangulation. The C i denotes the normalization constant when summing over index j ensuring that M is a right stochastic matrix [72].
We remark in our methods the approach of using the induced metric from the embedding allows for capturing the geometric contributions while avoiding the potentially cumbersome expressions that can arise from a more explicit treatment using differential geometry. Our methods utilize that the distances between points on the surface correspond to the arc-length determined by the path as measured using the embedding space and its associated notion of distance. As a result, in our methods the discretized surface inherits its metric without the need for further analytic derivations. Further properties of our methods include that the discretizations for curved surfaces are based on Markov-Chain transitions and as a consequence will have mass conservation up to numerical round-off errors.
The kernel M ij has been shown in the limit of refining the surface sampling to approximate diffusion under the surface Laplace-Beltrami operator D∆ M [81,82]. This has been shown to have the accuracy The O holds with → 0. The N is the number of points sampling M subject to a uniformity condition [42,81,82]. The u j = u(x j ) samples a smooth test function u. Intuitively, this follows since the stochastic matrix M converges to the operator as exp ε 2 ∆ M /4 = exp(D∆t∆ M ) exp(D∆t∆) and exp(D∆t∆ M ) I + D∆t∆ M I + D∆t∆, where ∆ = ∇ · ∇ is the standard Laplacian. The last two terms are motivated by the Taylor expansions of the exponential and the geometric terms as → 0.
Since the the surface properties will only be approximated when there are a sufficient number of sample points in the support of the kernel, for a given N there is a trade-off in the choice of . If is too large the approximation will not be of local surface properties. If is too small only the center point will contribute significantly to the kernel. We investigate further this trade-off in and the resulting approximation accuracy in Section II E.
We remark that our approach can also readily be modified to obtain models for proteins in super-sub diffusive regimes. One way this can be accomplished is by constructing random walks that are non-local on the mesh Shown is the relative error of the Markov-Chain surface discretization compared with the Laplace-Beltrami operator on the unit sphere S 2 . The case without area correction in equation (10) are shown as solid lines and with the area correction with dotted lines. There is a trade-off in to make it sufficiently large to contain enough points sampling the surface while also maintaining locality. For the truncation radius r0 ≥ 1.5 the optimum is around 1.
having long-tails or making state transitions between additional non-spatial states. These models could be constructed from theoretical considerations or empirical observations of autocorrelation statistics [76].

C. Detailed Balance and Area Corrections
To incorporate the contributions of the drift arising from U in equation (1), we consider the Gibbs-Boltzmann distribution expressed as ρ(x) = Z −1 exp (−βU (x)), with β = (k B T ) −1 the inverse thermal energy. We seek discretizations preserving statistical structure, such as detailed-balance as in [83,84]. For our curved surfaces, we seek discretizations for methods that have at steadystate the surface Gibbs-Boltzmann distribution with detailed balance [74]. We can express the evolution of the discrete probability in terms of net fluxes as At steady-state p k we design our transition rates so that the we have a discrete surface Gibbs-Boltzmann distribution with approximate detailed balance. The discrete detailed-balance J ij = 0 gives the conditions Motivated by these conditions, we discretize using the transition rates where C i normalizes the M ij to be a probability when summing over the index j. The conditions hold up to the normalization ratio C i /C j → 1 as the discretization is refined with N → ∞, → 0.

Algorithm 1 Hybrid Reaction-Diffusion Method
Update particles using jump rates from Section II C.
. Check for chemical reactions C and update the states.

D. Chemical Reactions and Field Coupling
We also formulate methods for capturing chemical reactions. We introduce for each particle X i (t) an associated "state" variable S i (t). The S i (t) can be thought of tracking the chemical species to which X i belongs at time t. We consider a collection of chemical reactions . For c we model the reaction using a Smoluchowski reaction radius r and probability p . For a second-order reaction, this means if two particles X n i and X n j have X n i − X n j < r , then we take probability p for a reaction to occur. For first-order reactions, we just take p = p (∆t) for each particle for the probability a reaction occurs over the time-step.
We also allow for coupling of the particles to underlying fields q = q(x, t). We discretize q in time and space using the lattice with q n a ∼ q(x a , t n ). As we discuss in later sections, the fields can influence the local jump rates and reactions. The evolution of the field is modeled using q n+1 = Ψ(q n ). The case of multiple fields can be handled by taking q to be vector-valued. We shall discuss in more detail specific evolution models for concentration fields in Section III A and phase-fields in Section III B. We give a summary of the steps in our method in Algorithm 1.

E. Approximation Errors
To obtain accurate diffusive dynamics we consider how parameter choices influence the jump rates and the approximation of the Laplace-Beltrami operator discussed in Section II C. We investigate the accuracy of the surface discretizations and the trade-offs in the choice of between sufficient sampling and maintaining locality. We perform our studies for the surface of the unit sphere S 2 , which has the Laplace-Beltrami operator ∆ M with eigenfunctions corresponding to the spherical harmonics [80,85,86]. In the comparisons, our numerical operator M is obtained from equation (11) and the scalings indicated in (8) to yield the approximation L = (M −I) ( 2 /4)∆ M . In practice for efficient calculations, the M ij is constructed by truncating the kernel only to use neighbors x j within the distance r 0 with |x i − x j | < r 0 . We consider the errors for a test function v given by In Figure 2, we show the relative errors based on rel = e where · 1 is the L 1 -norm averaging over the surface. We consider case when v is the spherical harmonic corresponding to v = v(x, y, z) = z restricted to the surface. The sphere is discretized using a triangular mesh with nearly uniform elements having N = 10, 0000 nodal points. We consider how the error varies for different choices of ∈ [0.25, 2.0] and r 0 ∈ [0.5, 2.0]. Letting δx = min ij |x i − x j |, we find that when δx there are insufficient number of points in the support of the kernel to estimate the surface geometry and the error becomes large. We find when δx is large there are many points within the support of the kernel, but the area of support is not localized enough to provide a good estimate of the operator. We also show both the case with and without the area correction terms. We find for our relatively uniform triangulations these give comparable overall errors here. For our discretizations of the sphere based on N = 10, 000 points and r 0 ≥ 1.5, we find that the optimal choice is 1, see Figure 2.

F. Role of Geometry in First-Passage Times and Other Statistics
We perform analysis to develop some results showing how our methods can be used for investigating the role of geometry of the first-passage times and other statistics associated with the drift-diffusion dynamics of particles on curved surfaces. Our Markov-Chain discretizations allow in some cases for computing efficiently statistics without the need to resort to Monte-Carlo sampling, provided the state space is not too large. We consider statistics of the form For the domain Ω, the τ Ω = inf{k ≥ 0 | X (k) ∈ Ω} is the stopping time index for the process to reach the boundary ∂Ω. Each of these statistics can be computed without sampling by the following results.
Proof. (see Appendix A) Theorem 2. The statistics w (Ω) in equation 14 satisfies The ∂w extracts entries for all indices with x i ∈ ∂Ω. Thê M ,Î refers to the matrix only with the rows with indices in the interior of Ω. The surfaces have area one and when there are edges we use reflecting boundary conditions. For the shapes (i)-(ii) spotted patterns emerge having roughly a hexagonal pattern. For spherical topology (iii) a regular hexagonal pattern without defects is no longer possible, and instead striped patterns mix with spots. For the case of a torus (iv), which can sustain a hexagonal pattern in principle, the heterogeneity of the curvature appears to drive the formation of localized stripe-like patterns. These results indicate that both the geometry and topology can significantly impact pattern formation.

Proof. (see Appendix A)
In the case that f = 0, g = 1, this becomes the First- These results allow for the statistics of equation 13 and 14 to be computed efficiently without the need for Monte-Carlo sampling provided the state space is not too large.
Motivated by observations that the neck geometry appears to be a strong factor in compartmentalization in dendritic spines [68,87], we compute the first passage times of protein diffusions for different spine shapes in Figure 3. The geometries were generated from isosurfaces in a technique commonly referred to as meatballs [88]. We used several spheres to form the tubular region and another sphere for the bulb region. The geometry was obtained as the isosurface using the level set function of the form f (x) = k φ k (|x − x k |), where φ(r) = c k /r 2 . The surface was triangulated at a refined spatial resolution to obtain around N = 29, 000 nodal points for each shape.
First passage times are computed for when a protein starts at the top of the head region with the bulb-like shape and reaches the boundary. To investigate the role played by the neck region compared to the influence of the other aspects of the geometry, three cases are considered for the boundary location. These are boundaries located (i) at the base of the spine, (ii) in the middle of the neck, and (iii) within the spherical head just above the neck. These results indicate that compared to the other geometric features, the width of the neck region plays the dominant role in the first passage times. As the neck narrows, the first-passage time for the protein diffusion significantly increases from these changes in the geometry, see Figure 3. Shown is the evolution of the reaction-diffusion pattern formation process on the dendritic spine shape having the narrowest neck, label A (notation the same as in Figure 3). The pattern progresses first by forming within the bulb-like head region, and then spreads through the neck to form patterning on the tubular part of the domain. The time-steps are shown for n = 0, 400, 800, 1200, 1600, 2000.

A. Turing Instabilities on Curved Surfaces
We show how our discretization approach can be used to develop methods for performing simulations of general reaction-diffusion processes on curved surfaces of different shapes. We consider reaction-diffusions, such as the pattern formation process based on Turing's instability mechanism [89], where the geometry and topology of the domain can impact the patterns that are obtained [90,91]. Consider the system with two molecular species with concentrations u, v The diffusivities D u and D v will in general be different. Through the non-linear reaction terms, the difference in diffusivity can cause the homogeneously mixed concentrations to become unstable resulting in pattern generation [89,90]. We consider Gray-Scott reactions [92], which can exhibit different patterns depending on the initial conditions and interactions with noise and other perturba- tions [90,93]. For the Gray-Scott reactions [92], the terms are f (u, v) = −uv 2 + a(1 − u) and g(u, v) = uv 2 − (a + b)v. The rate parameters a, b are for the chemical reactions u + 2v a − → 3v and v b − → p. We start with a homogeneous steady-state solution for the system (u, v) = (1, 0), and add small perturbations based on uniform random noise ±10% of the steady-state at each location x i within a region Γ ⊆ Ω. Throughout our simulations, we choose a = 0.04, b = 0.06, guided by the parameters of the phase diagram of [94]. We investigate how the Gray-Scott patterns are influenced by different geometries and topologies.
Using our Markov-Chain discretizations, we develop numerical methods for the spatial-temporal evolution of the concentration fields. By equation (7), we obtain a stochastic matrix M . The overall idea is to model the concentrations by u(x, t) = (n u /A)p u (x, t), v(x, t) = (n v /A)p v (x, t), where A is the surface area. The p = [p u , p v ], with [p n ] i = p(x i , t n ), are obtained from the probability evolution for the Markov-Chain given by p n = p n+1 M . In this way, we obtain a model that approximates the diffusive evolution of the continuous concentration fields.
To model the full reaction-diffusion system in our simulations, we split the time-step integration into a diffusive step and a reaction step. For the diffusive step ∆t, we use our Markov-Chain discretization in Section II C to update the concentration fields q = (u, v). This provides a mean-field model for the concentration field, with q n+1/2 = Ψ 1 (q n ) corresponding to u n+1/2 = u n M, v n+1/2 = v n M . The M is the transition matrix for the surface diffusion derived in Section II C. For the reactions, we use the smaller time-steps δt = ∆t/n rk which are integrated with the fourth-order Runge-Kutta method RK4 [95,96]. This gives q n+1 = Ψ 2 (q n+1/2 ). We alternate between these steps in our simulations to obtain q n+1 = Ψ(q n ) = Ψ 2 (Ψ 1 (q n )).
As we discussed in Section II B, the M approximates the Laplace-Beltrami operator. This ensures with enough spatial resolution the evolution will approximate diffusion on a surface. One way to view our model is as a meshfree way to discretize the surface reaction-diffusion PDEs. This provides simulation methods corresponding to the mean-field concentration fields of the protein species.
Our approaches for the diffusion also can be used more generally to go beyond the mean-field model by using discrete particle simulations with individual random walkers. These would have the same distribution as our continuum model when taking appropriate limits. Our methods allow for either (i) to perform stochastic simulations of random walks tracking individual particles to account for discrete spatial-temporal fluctuations from finite number effects, or (ii) to perform deterministic simulations tracking the probability distribution of the walkers. For the reaction-diffusion studies we use approach (ii). For later studies involving phase separation, we use approach (i).
We consider the Gray-Scott reaction-diffusion process on the following geometries (i) flat sheet with boundaries, (ii) finite cylinder with boundaries, (iii) surface of a sphere, (iv) surface of a torus, see Figure 4. The surfaces have area one and when there are boundary edges we use reflecting boundary conditions. For the shapes (i)-(ii) spotted patterns emerge having roughly a hexagonal pattern. For the spherical topology (iii) a regular hexagonal pattern without defects is no longer possible, and instead striped patterns mix with spots. For the case of a torus (iv), which can sustain a hexagonal pattern in principle, the heterogeneity of the curvature appears to drive the formation of localized stripe-like patterns. The results obtained with our methods indicate that both the geometry (curvature and scale effects) and the topology can impact significantly the pattern formation process.
We remark that Turing instabilities as a mechanism for pattern formation were originally motivated by a linear stability analysis on periodic domain performed by A. M. Turing in [89]. Extending the analysis to general manifolds is non-trivial given the challenges of analyzing even for the linearized dynamics the eigenfunctions and eigenvalues associated with the Laplace-Beltrami based-diffusions and roles of curvature and topology. Some work in this direction has been done in the literature in the case of specialized geometry, such as the sphere [97,98]. The results of our methods for Turing instabilities yield patterns comparable to those seen in such previous studies. In the spherical case we find that there is a transition from striped to spotted patterns in the sphere size [99]. We also find that the hexagonal spotted patterns with defects manifest with arrangements similar to [97,98].
We also consider more complicated shapes for the Gray-Scott reactions given by our mechanistic dendritic spine geometries with different neck sizes, see Figures 5,6. These shapes consist of a bulb-like head region (representing the spine) which is connected by a neck structure to a tube-like region (representing the dendrite). We vary the shapes by changing the thickness of the necklike structure joining the two regions. The evolution of the pattern formation process for the geometry with the narrowest and widest neck shapes are shown in Figure 5. For the narrowest necks, the confinement in the head region appears to result in stripe-like patterns that also extend through the neck. In the larger tube-like region, the spotted patterns form inter-mixed with the stripe pattern, see Figure 6. These results indicate how local regions can exhibit different patterning depending on the local geometry.

B. Dendritic Spines and Protein Kinetics
We develop a mechanistic model for dendritic spines to investigate the dependence of protein transport and kinetics on geometry and heterogeneities arising from phase separation. Our investigations are motivated from exper-  iments on SynGAP and PSD-95 proteins, where phase separation may play a role in driving receptor organization [70]. Phase separation can arise from nucleation or modulating local concentrations in the cell membrane, which is a heterogeneous mixture of lipids, proteins, and other small molecules [1,100,101]. In our mechanistic model, we investigate how the spine geometry and phase separation can influence reaction kinetics. We start with a two species system. The A and B molecular species are tracked at the individual particle level. In the absence of phase separation, this would diffuse freely over the membrane surface. As a starting point, we study the basic chemical kinetics A + B → C. The discrete particles react with probability p when coming within a reaction distance r < r 0 . This is motivated by Smoluchowski reaction kinetics [102][103][104][105].
To model such effects at a coarse-level, we track a continuum phase field φ = φ(x, t) which can couple to the local diffusive motions of the discrete protein particles. We consider in our initial model the case of a local order parameter associated with phase separation based on Ginzburg-Landau (GL) theory [106,107]. Other phaseseparation phenomena and approaches also could be considered in principle within our framework, such as the second-order Allen-Cahn [108] or the conservative fourthorder Cahn-Hilliard [109], with additional work on the numerical methods for the operators on curved surfaces and coupling with the Markov-Chain discretization [110]. For notational convenience, we will denote q φ(x), so q : M → R for the phase variable map for the curved surface M. This gives at lattice site i, The first term accounts for a line tension between phases and the second term drives the phase to ±1 with local energy To capture similar effects aŝ V 0 [q], we use the simplified discrete model with The first term models the interfacial tension, where the coefficient The second term models the order parameter which locally is at a minimum for q = ±1. The α 1 , α 2 > 0 control the strength of the interfacial tension and the phase variable ordering. We use for the decay coefficient W ij = exp −r 2 ij /ε 2 which is truncated for large r ij . To couple protein motions X(t) and their ability to induce local phase ordering in q, we use the free energy term The α 3 > 0 controls the strength of this coupling. The η a gives a kernel function with support localized around the protein location X = X(t). This term drives the phase field toward q = −1 in a region around the protein location X(t). We use η a (r) = exp −r 2 /a 2 and The free energy for full system (q, X) with protein configuration X and phase-field q is The discrete protein positions [X] k = X k are updated using the Markov-Chain discretization in equation (11) FIG. 9. Dendritic Spine Model: Average Phase. The proteins C can nucleate phase with q = −1. Shown are how the geometry impacts phase evolution for five different sizes for the neck region with labels A-F ordered from narrowest to widest as in Figure 6. The phase separation is found to slow down considerably for the narrowest neck (label A) relative to the widest neck (label F). The differences emerge around around time-step ∼ 300, when the bulb-like head region is saturated. We see over the 2000 time-steps the widest neck results in almost the entire domain including the tubular region converting to phase q = −1. For the narrowest neck shape, we see at time-step 2000 the phase separation is primarily isolated to the head head region and surrounding area, while the rest of the domain remains primarily q = 1. The average is taken over the entire surface.
for the drift-diffusion dynamics. The energy for the protein configuration is given by U (X) = V [q, X]. The time evolution of the phase field is given by This is discretized and integrated over sub-time-steps δt = ∆t/n rk using Runge-Kutta RK4 [95,96]. Modeling the discrete system using the common free energy V ensures the bi-directional coupling gives equal-and-opposite forces between the phase-field q and the proteins X.
We perform simulations to investigate how the dendritic spine geometry impacts the protein reaction kinetics and phase-separation, see Figure 7 and Table I. The proteins couple to the phase separation by having the ability to drive nucleation of local patches with q = −1 near the protein location X(t). The geometries were chosen with the neck region taking on shapes varying from narrow to wide, see Figure 4. The protein species A, B are modeled as originating in the top head region of the spine. This could arise for instance if these species correspond to tracking proteins only after they have become activated in this region. The proteins can then diffuse, and may induce local phase-separation through the coupling V 3 . The proteins can also interact to form a complex which is tracked by molecular species C.
We show how the phase separation proceeds over time as the size of the neck region varies in Figure 8. As the neck region narrows, the protein diffusion and the phase separation are restricted to be localized in the head region. As the neck becomes wide, the protein diffusion and phase separation can readily proceed to spread more rapidly into the tubular region, see Figures 8,9. Further effects arise from the coupling of the proteins with the local phase.
The protein dynamics are impacted by both the geometry and the local phase. In our simulation studies, the coupling coefficients are taken to be α A 3 = α B 3 = 0 and α C 3 = 1.0. For this case C can nucleate phase separation nearby since it prefers the phase field have q = −1. This also results for C experiencing a trapping force within regions with q = −1, from the free energy in equation (22). This restricts the movements of C. The evolution of the creation of of molecular species C is shown in Figure 10. The number of C formed and retained in the head region is significantly impacted by the neck geometry and phase separation. When the neck is narrow, the confinement of A, B and the phase separation to the head region, both enhances the creation of C by more frequent A − B encounters and in C's retention to the head region from the phase trapping forces. As the neck becomes wide, the A, B can diffuse more freely and when the phase does nucleate it can more rapidly spread throughout the whole domain. This results in a much smaller number of C being created and retained in the head region.
Our basic mechanistic model shows that the morphology and phase separation can interact to serve together to enhance retaining proteins near the top of the dendritic spine. These results are expected to carry over to more complicated chemical kinetic systems. The further interactions between diffusion, kinetics, phase separation, and the geometry can regulate in different ways the spa- tial arrangements and the local effective rates of reactions in curved cell membranes.

C. Protein Clustering and Role of Phase Separation
We investigate how phase separation can impact protein kinetics and clustering. Consider the case of states where proteins in state A are in an active state receptive to binding and proteins in state B are dormant. Consider the competing reactions for forming dimers The p and q denote here the reaction probabilities. The protein-phase-field interactions are modeled using our approach in Section III B. We take the phase field to impact the proteins in states A and 2A with the same strength, In contrast, we take proteins in the dormant state B to not interact with the phase field, α B 3 = 0. To study the effects of the phase field, we vary α A 3 over the range from 0.0 (no interaction) to 1 × 10 1 (strong interaction). We investigate in simulations the number of dimers 2A proteins that are formed over time for different levels of phase field interaction α A 3 and as the phase separation progresses. The 2A counts are given in Figure 11. We find in the strongest coupling cases the phase field can significantly impact the reaction kinetics promoting clustering and dimer formation. The details of the parameters used in our simulations are given in Table IV. The results show that phase-field separation can play significant roles impacting reaction kinetics within membranes. In this case promoting the formation of clusters. The underlying mechanisms has to do with the formation for proteins of local environments with initially small phase-separated domains that serve to confine the diffusion of the proteins. Relative to free diffusion over the entire membrane surface, the phase confinement results in diffusion with the phase domain and more frequent encounters between co-confined proteins. In this way, the binding kinetics are enhances relative to the free diffu- sion. From the greater slope in the increase in the number of clusters formed for case α A 3 = 1 × 10 1 vs α A 3 = 0, we see the phases have the greatest impact when there are small phase domains. The smaller domains results in more frequent encounters between confined proteins. As the phase separation progresses the domains become larger and merge with the effect of the phase-coupling on the reactions become less pronounced.
These results show how phase separation and domain formation can be utilized within cell membranes to significantly augment the reaction kinetics without changing the reaction characteristics of the individual proteins. These simulation methods provide further ways to investigate the collective kinetics of proteins and the roles played by phase composition and the membrane geometry.

IV. CONCLUSIONS
We have developed biophysical models for investigating protein kinetics within heterogeneous cell membranes of arbitrary shape. The introduced methods allow for studying at the continuum and discrete protein-level the drift-diffusion dynamics and reaction kinetics taking into account the roles of geometry and coupling to membrane phase separation.
We have shown that coupling of kinetics to phase separation and geometry can significantly enhance or inhibit reactions and cluster formation processes. These initial results suggest a few mechanisms by which biological membranes may utilize phase separation and geometry to drive kinetics and organization of proteins within cell membranes.
We discretize the particle drift-diffusion dynamics on the surface using a Markov-Chain with jump rates M ij . The First Passage Time (FPT) and other statistics can be computed efficiently from the Markov-Chain without the need for Monte-Carlo sampling when the discrete state space is not too large. We show how results similar to the Backward-Kolmogorov PDEs [75] can be obtained for our discrete Markov-Chains. Proof. For the initial condition X (n) = x i , let the matrix P (m) have the components P where τ Ω = inf{k ≥ 0 | X (k) / ∈ Ω}. The w satisfies the linear equation The ∂w extracts entries for all indices with x i ∈ ∂Ω. Thê M refers to the rows with indices in the interior of Ω.
Proof. Since the equation is linear, we will first consider the case with f = 0, g = 0 and then the case with f = 0, g = 0. The general solution is then the sum of these two cases. For the case with f = 0, g = 0 we have g(X (t) ) X (0) = x i , τ = n Pr{τ = n} TheÎ is the linear map that extracts entries within the interior of the domain Ω.
We next consider the case with f = 0 and g = 0. Similarly, this follows from For x i ∈ ∂Ω we have w i = E f (X (τ ) ) X (0) = x i = f (x i ), since τ = 0 in this case. In matrix form this gives ∂w = f and (M −Î)w = 0. (A12) The ∂w extracts the entries with indices i corresponding to the boundary ∂Ω. Putting both cases together we have that w satisfies for general f , g the linear system These results provide approaches for computing efficiently the statistics u and w without the need in some cases for Monte-Carlo sampling when the state space is not too large. These results provide an analogue of the Backward-Kolmogorov PDEs in the setting of our Markov-Chain discretizations of particle drift-diffusion dynamics on curved surfaces.