Physical phase field model for phagocytosis

We propose and study a simple, physical model for phagocytosis, i.e. the active, actin-mediated uptake of micron-sized particles by biological cells. The cell is described by the phase field method and the driving mechanisms of uptake are actin ratcheting, modeled by a dynamic vector field, as well as cell-particle adhesion due to receptor-ligand binding. We first test the modeling framework for the symmetric situation of a spherical cell engulfing a fixed spherical particle. We then exemplify its versatility by studying various asymmetric situations like different particle shapes and orientations, as well as the simultaneous uptake of two particles. In addition, we perform a perturbation theory of a slightly modified model version in the symmetric setting, allowing to derive a reduced model, shedding light on the effective driving forces and being easier to solve. This work is meant as a first step in describing phagocytosis and we discuss several effects that are amenable to future modeling within the same framework.


Introduction
Biological cells often have to engulf larger particles [1], either because they feed on them or as part of the immune response against intruders like bacteria, or dying cells.While the cell membrane's primary function is that of a barrier, it has to be much more versatile and dynamic [2] and hence cells have developed several mechanisms to ingest particles: uptake of small (up to ∼ 100 nm) particles is called endocytosis and is typically passive, mediated by receptor-ligand binding or clathrin self-assembly [3].In contrast, larger micron-sized particles are actively taken up by a process called phagocytosis.Finally, there is a third and less specific way, called macropinocytosis, where membrane protrusions are actively formed "out of the blue" and the cell "swallows" fluid that may contain several particles [4].
Phagocytosis modeling has been reviewed recently in Ref. [5].Models mostly focus on the so-called zipper mechanism [6] of receptor-ligand binding, overcoming membrane restoring forces.In fact, in most phagocytes, Fc receptors (FcR) in the cell membrane are attaching to immunoglobulin G (IgG) present on the particle's surface, leading to a zipper-like advancement of a thin membrane protrusion around the particle.Refinements of such models also account for signaling and coupling to active processes, phenomenologically modeling the effects of actin.On the other hand, a framework focusing more on flows and forces, and employing a two-fluid model, for the cytoplasm and the actin network, was developed in Refs.[7,8].Very recently, a Monte-Carlo approach for the membrane shape, including "active energy" from actin, was proposed as well [9].
While it was clear early on that self-assembling, biochemically driven actin filaments play some role [10], recent experiments showed more explicitly than ever that particle internalization in phagocytosis is driven by the formation of actin protrusions that wrap around the particle [11].It is therefore crucial to develop models including the actin dynamics underneath the deformable and dynamic membrane close to the particle.Upon closer inspection, it turns out that phagocytosis and cell motility are quite related to each other: in fact, both rely on actin-driven membrane protrusions, on the adhesion of these to a substrate (round particle vs. often flat substrate, FcR-IgG vs. integrin binding to fibronectin or related ligands), and confront us with the numerical problem of moving and deformable boundaries.Luckily, the modeling of cell motility has seen much progress in recent years, mostly due to the phase field framework [12].Just recently, a phase field approach has already been applied to macropinocytosis [13], however implementing activator-inhibitor fields instead of actin.
We here propose a simple, physical modeling framework to model the initial stages of phagocytosis.It is based on a three-dimensional (3D) phase field model for crawling cells [14], only slightly adapted to the problem under consideration.We show the framework's versatility by studying the uptake dynamics of spheres as a function of actin ratcheting, the effects of particle shape and orientation, and the simultaneous uptake of two particles.In addition, we use perturbation theory to derive a reduced equation for the cell shape, valid for the initial regime of phagocytosis, that is easy to interpret and faster to solve.

Basic model framework
We base ourselves on a phase field modeling framework that was originally developed for cell motility, first in 2D [15] and then in 3D [14].We stick here to the basic version of this framework, to make the model as simple and transparent as possible.Importantly, however, many model extensions and variations have already been included in the context of cell motility, for instance, adhesion formation to inhomogeneous substrates [16], flows [17] and elastic effects [18] in the cytoskeleton, as well as regulation pathways [19] and activator-inhibition dynamics [20], which all pave routes to future refined models for phagocytosis as well.
The approach is conceptually simple and flexible and has two dynamic variables: the cell's shape is described by a dynamic 3D scalar phase field ρ(r, t) ∈ [0, 1] and its transition region is identified with the cell membrane (ρ = 1 corresponding to the cell and ρ = 0 to the outside).The second dynamic variable is a 3D vector field p(r, t) that "lives" only inside the domain defined by ρ and that describes the mean degree of parallel ordering (absolute value) and the mean orientation (vectorial direction) of the structurally and dynamically polar actin filaments in the cytosol.Additional phase field(s), that can be either static or dynamic, can then be used to define regions with steric exclusion, allowing to implement arbitrary substrates or environments [14].In the case of interest here, we will use such an additional phase field Φ(r) to describe a particle of a chosen shape that may get close to the cell and internalized via phagocytosis.
The actin dynamics inside the cell is implemented on the phenomenological level via its basic physical features: actin is nucleated close to the membrane (i.e. the phase field interface) and exerts a force against the membrane via the well-known polymerization ratchet mechanism [21], pushing it forwards.Two more properties of actin are important: i) its nucleation is localized also close to the substrate/particle since its regulators receive signals there.This is implemented via another phase field Ψ(r) (which, if not stated otherwise, we choose here to be the same as Φ(r) for simplicity) having a certain decay length away from the substrate/particle ‡ and as such defining a region where actin nucleation is possible [22].And ii), actin is oriented predominantly tangentially to the local substrate/particle surface, which can be modeled via a projection operator onto the local tangential plane.
Explictly, the 3D model equations read with the two static fields Φ(r), Ψ(r) defining the particle, and the localization of its stimulation of actin polymerization close to that particle, respectively.The first two terms in Eq. ( 1) describe classical phase field dynamics [15,12]: D ρ fixes the width of the diffuse interface and the second term implements a relaxational dynamics in the phase field double well potential.Since in this simplest implementation the dynamics of ρ is non-conserved, it needs to be supplemented by a volume conservation, implemented here via a global constraint, δ , with the second term penalizing differences between the cell's current volume ρ d 3 r and the prescribed volume V 0 , the stronger the larger the "stiffness" µ.
The last three terms in Eq. ( 1) model processes that can move the phase field boundary: αp • ∇ρ implements the pushing by actin ratcheting against the cell's boundary, with α an effective velocity of the interface movement [15].κ∇Φ • ∇ρ ‡ Both Φ(r) and Ψ(r) are, in the simplest case of a spherical bead, implemented as a function with a tanh-like radial profile around the particle's center, see also Appendix A for more information about the implementation.
implements adhesion [23,24] between the cell and the particle described by Φ(r) by means of ligand-receptor bond formation, with κ the adhesion strength parameter.Finally, the last term models excluded volume interaction with the particle; it can be derived from an interaction potential λ 2 ρ 2 Φ 2 d 3 r with excluded volume strength λ.For more information on these contributions see Ref. [25].
The actin dynamics, Eq. ( 2), has three contributions from actin turnover: a diffusion/elastic term, a source term proportional to the polymerization rate β and a sink term proportional to τ −1 describing depolymerization, with τ of the order of the actin turnover time.Since the actin creation and its pushing via the polymerization ratchet are closely related, we typically fix β = 3α, which was a reasonable value when applying the model in the context of cell motility [14].The geometric structure of the actin source term was developed in Ref. [14]: actin polymerization is modeled anisotropically with respect to the local tangential plane defined by Φ(r) modeling the particle.This is achieved by employing the projection operator onto that plane, P = Î − n⊗ n with Î the identity and n = ∇Φ/|∇Φ| the local normal vector on the particle's surface.Finally, the last term in Eq. ( 2) is added to ensure that there is no actin penetrating into the particle, which would be unphysical.
It should be noted that the used phase field approach has an inherent surface tension associated with the phase field's wall energy [12].This apparent surface tension is related to the excess energy due to the phase field interface and is Σ ∝ D ρ .The surface tension could be removed [26] or tuned to a desired value [27], but we keep it here, since it reflects the only restoring forces due to membrane shape changes (since we do not consider membrane bending yet), which is a relevant effect in phagocytosis.
Details on the numerical implementation of the model can be found in Appendix A with typical parameter values given in table A1.Fig. 1 shows applications to several example geometries relevant for phagocytosis, highlighting the versatility of the modeling framework: A1 and A2 show the engulfment of a spherical particle of small size, while in B1 and B2 the particle is of similar size as the phagocytosing cell.C shows the initial phase of the uptake of a spherocylindrical particle.D exemplifies the simultaneous uptake of two small spheres and E the uptake of a sphere by a cell that is spreading on a flat substrate.

Results for a spherical cell engulfing a spherical particle
We first applied the model to the most symmetric and most studied situation: an initially perfectly round, freely floating cell that engulfs a static, non-movable, round bead, cf.Fig. 1 A,B.This situation has been experimentally studied, e.g., in Refs.[7,8] using a doube-micropipette apparatus and is often used as a benchmark for more complex geometries and scenarios [5].
To quantify the uptake dynamics, we introduce the following observables: first, the relative cell-particle contact area, A(t) ∈ [0, 1], which we define in the phase-field sense as the integrated overlap of the cell phase field, ρ, and the static phase field describing Here, the first phase with the formation of a phagocytic cup is clearly visible, followed by a slower relaxation, where the whole cell body rounds up.Snapshots correspond to t = 50s (A1, B1) and t = 400s (A2, B2).C) A spherocylinder (radius r = 4; cylindrical element of length 4r) is engulfed in a "tip-first" configuration; t = 15s.Sufficiently high actin activity leads to full engulfment after about t = 400s.D) Simultaneous engulfment of two particles (of size r p = 4), positioned at an angle of 90 • with respect to the initial center of the cell; early stage at t = 10s.E) A cell was allowed to spread on a flat substrate (yellow) modeled by an additional static phase field.Then a small sphere was approached.Parameters: A, B, D: α = 4; C: α = 6; others as given in table A1.E: as specified in the main text, see section 4.
the particle, Φ: Here A f ull is the overlap for full engulfment, obtained by a reference simulation placing a particle in the center of the cell for the same parameters.Clearly, A ∈ [0, 1] holds, with A = 1 meaning full engulfment.
To also quantify the mechanism of phagocytosis, as a second observable we measure the total amount of actin generated in the cell, normalized by the total cell volume: Note that actin is only nucleated in a region close to the particle §, defined by the field § In the case of a cell that in addition spreads by actin-polymerisation, cf.Fig. 1E, one would need to B A C Figure 2: Quantification of the uptake dynamics of spherical beads by initially spherical cells.A) shows the relative cell-particle contact area, A(t), for r c = 14 and r p = 10 and for different actin pushing rates α.For sufficiently large α, there is a fast initial dynamics (black line indicates t 2 ), followed by a second slower regime.B) shows the total amount of actin, P (t), for the same situation as in A. C) Shows the final engulfment A(t → ∞) color-coded as a function of actin strength and particle radius.With increasing actin polymerization and ratcheting parameter α, particle uptake becomes faster and more complete.The size dependence is weak, since larger particles also stimulate more actin growth.Parameters as given in table A1.
Ψ(r), hence P (t) directly quantifies the amount of actin in the phagocytic cup.
Uptake dynamics.Fig. 2A shows the relative cell-particle contact area as a function of time, A(t), for a cell of radius r c = 14 engulfing a large particle of r p = 10 (cf.Fig. 1 panels B1, B2 for two snapshots) and varying actin pushing rate α.Note that the actin polymerization rate then also increases since we chose β = 3α.For sufficiently large α, one can see an initial fast dynamics (the black line shows the power law A(t) ∝ t 2 as a guide to the eye) associated with the cup phase.This is followed by a second regime of much slower engulfment dynamics.For smaller α the uptake is slower overall, and the two regimes are less distinct.Engulfment dynamics has been carefully analyzed in Ref. [28], using data from Ref. [8], and interpreted using a receptor-based model.There it was found that the engulfed arc length scales like √ t for diffusive receptors and linearly in t for perfect receptor drift.Our model reflects the second case -note that we measure the area of contact, A(t) ∝ t 2 , which scales like arc length squared -due to the constant pushing of actin.Experimentally, the behavior is much more complex, with an apparant jump between a diffusive and an imperfect drift regime, that could be modeled only by using a signaling molecule in addition [28].
Fig. 2B shows the volume-normalized total amount of actin, P (t), that is generated in the cell due to the interaction with the particle.Note that only in a region close to the particle (modeled via the static phase field Ψ) actin generation is stimulated.Its dynamics is characterized by the actin-related parameters α (reflecting ratcheting) generalize the definition of P (t) accordingly.and β (nucleation).After the initial cup phase, the actin forms a ring-like structure at the protrusion front that moves forward and engulfs the particle.For sufficiently large α, P (t) hence increases in time, reaches a maximum at roughly half engulfment, and then decreases again.This behavior is absent for small α values.The maximum in P (t) is mostly a geometric effect, since the actin ring forming around the particle at the advancing front is largest at half-wrapping.
Fig. 2C shows an "uptake diagram" as a function of actin ratcheting strength α and particle radius r p (for fixed cell radius r c = 14).The color code marks the "final engulfment", A(t → ∞), obtained by running long simulations and extrapolating.It is clearly visible that increasing α fosters engulfment, cf. also supplementary movies 1 & 2. The particle size dependence is rather weak, with larger particles being engulfed slightly easier.This is due to the fact that the particle stimulates the generation of actin, hence a larger bead directly translates to more actin (cf.also Fig. 2B).The simple model proposed here so far neglects effects like limited resources for actin or time-delay due to actin having to diffuse into the engulfing zone, which could be included in a future, more refined model.
We also investigated the effect of surface tension.As discussed in section 2, the phase field method has an inherent surface tension associated with the phase field's wall energy.It can be suppressed by using an additional term like +D ρ c|∇ρ|, with c = −∇ • (∇ρ/|∇ρ|) the curvature of the cell's phase field interface, in Eq. ( 1) [26,27].Doing so, we generally observe an enhanced engulfment capability of the cell (data not shown).The effect is most pronounced in an intermediate parameter regime, where shape changes are already significant but the active forces are not large enough to quickly and fully engulf the particle.This confirms that, as expected, surface tension comprises a restoring force to interface deformations and hence counteracts phagocytosis.

Select asymmetric problems
Effect of particle shape.Phagocytosis is known to depend on particle shape [29,30,31].While spherocylinders, cf.Fig. 1C, are relevant since bacteria often have these shapes, the most studied and instructive shape variations are rotational ellipsoids (spheroids).Experimentally, in Refs.[32,33] it was found that oblates are engulfed easier than spheres, which in turn are easier to engulf than prolates.In addition, for non-spherical particles the orientation in which it meets the cell also is important.In [29] the simple idea was put forward that the curvature at the initial point of contact determines -at least the initial phase of -phagocytosis.
Fig. 3 compares the uptake dynamics of a sphere, an oblate ellipsoid ("lentil") and a prolate ellipsoid ("rugby-ball"), with the ellipsoids in the two possible main orientations.The volume of all particles was chosen equal, and the aspect ratios of the elliposids was chosen to be 1/4 and 4, respectively.In Fig. 3A we consider the case where receptorligand mediated adhesion (κ = 12) dominates over actin ratcheting (α = 1).One can see from the dynamics of the relative cell-particle contact area A(t) that less curvature

A B C
Figure 3: Effect of particle shape on uptake dynamics.A) shows the relative cellparticle contact area, A(t), as a function of time for a sphere (black), for an oblate ellipsoid touching the cell with its bottom (violet) and with its rim (blue), and for a prolate ellipsoid touching the cell with its side (red) or tip (brown).All particles have equal volume.Parameters are such that there is strong receptor-ligand mediated adhesion (κ = 12) but weak actin ratcheting (α = 1).B) Same as A, but now the respective A(t) of the sphere is substracted.C) Same as B, but now for larger actin ratcheting (α = 4).Bottom: Snapshots for the respective shapes and orientations.Shown is the case with α = 4 at t = 10.
is beneficial for the initial (so-called pedestal) phase.The uptake processes are in the following order: oblate (bottom) has the lowest curvature and is fastest, followed by the prolate (side), the sphere, the oblate (rim) and finally the prolate (tip) which has highest curvature and is slowest.This is in accordance with receptor-based models, which show that higher-curved regions are harder to engulf [5].For better comparison, Fig. 3B also shows the deviations from the A(t) of the sphere.Fig. 3C shows the deviations from the A(t) of the sphere for the case where actin ratcheting is more important (α = 4).Interestingly, the initial phase has the same ordering, but then the least curved oblate (bottom) slows down substantially and the prolate (side) slightly.In turn, the oblate (rim) now performs much better, as does the prolate (tip).This clearly indicates that the adhesion dynamics and the actin-ratcheting dynamics have different preferences concerning curvature.Finally, the bottom row of Fig. 3 shows snapshots of the initial stage of engulfment for the various shapes and orientations, cf. also supplementary movies 3 & 4.

Simultaneous phagocytosis of two beads.
It is known that cells can engulf more than one particle at the same time.For several reasons -among them, the increase in membrane tension, limited resources, as well as signaling -one would expect collective  A1.
effects to occur.Recently, a spatial resolution limit for phagocytosis was discussed in [34], using holographic optical tweezers to approach two IgG-coated beads with controlled distances to macrophage cells.While signaling is not included in our modeling framework yet (except for local actin stimulation ∝ β), we can explore the collective effects arising from geometry and volume conservation.Fig. 4 shows results for the simultaneous uptake of two spherical beads (r p = 4) by an initially spherical cell.The beads were arranged such that they touch the cell at an angle of 45 • , 90 • and 180 • with respect to the cell's center.Panel A shows the difference in contact area with respect to the uptake of a single sphere, ∆A(t) = A two (t)/2−A single (t).Interestingly, the 45 • positioning (green curve) leads -after an initial small slowdownto faster uptake than in the single bead case, before the difference decreases for longer times as A → 1.In contrast, the other two configurations lead to slower uptake and to less overall engulfment (blue and red).Panel B shows the actin observable P (t), measured per bead, showing the same trend: the 45 • configuration initially needs more time but than is faster than for single bead, while the other two settings get "stuck" at higher level of P , reflecting incomplete engulfment.Panels C1-C3 show snapshots for the three configurations at equal times, t = 50s.
Phagocytosis by cells spreading on a substrate.In dedicated experiments, phagocytosis can be probed using a freely floating cell sucked into one pipette while presenting the particle with another pipette [8].In practice, however, most experiments study cells that already spread on a substrate when phagocytosis is initiated.This means that the initial symmetry of the cell is already broken.We also studied this situation within our proposed framework, cf.Fig. 1E.The cell was allowed to spread on a substrate (yellow, modeled by another static phase field).Then a spherical particle was presented from above and engulfment started.
It should be noted that, in general, both receptor-ligand interaction and actin nucleation and polymerization dynamics are different for the cell-particle (κ, β) and cell-substrate (κ s , β s ) contacts.Due to the use of different phase fields for the particle and the substrate, this can be easily implemented: Fig. 1E shows the case of a cell (r c = 14, α = 2) that is moderately spreading on the substrate (κ s = 5, β s = 6) and has a strong interaction with the particle (κ = 18, β = 12) of size r p = 4.
We found (data not shown) that the uptake is slower and more incomplete compared to an initially spherical cell, since the spreading already increased the cell's surface area, leaving less freedom to engulf particles.We also studied the simultaneous uptake of two spherical beads, presented to the spreading cell in a plane parallel to the substrate, and found the same trends as discussed in the previous section, see supplementary movie 5, i.e. accelerated uptake for close-by beads (45 • ) and slowed down uptake for larger angles.

Reduced description of uptake dynamics from perturbation theory
One advantage of having a continuum PDE model at hand is that it is amenable to perturbation theory to derive reduced models.These typically have a restricted validity range, but are easier to solve numerically and yield additional analytical insight.For instance, from a full phase field framework, reduced equations for membrane waves [35,36], and cell spreading [37] have been derived in 2D and 3D, respectively.
The geometry considered is the one discussed in section 3, i.e. an initially perfectly spherical cell touching and starting to engulf a spherical bead, see Fig. 5A for a sketch of the geometry and definition of the variables.The approach is detailed in Appendix B. In brief, a perturbation method is applied in the sharp interface limit [38]: one assumes that ϵ, the ratio of the cell's phase field interface thickness to the cell radius, is small, which is typically well fulfilled.The model equations ( 1) and (2) had to be slightly modified concerning the definitions of Φ and Ψ, to avoid the perturbation to become singular.Assuming a suitable scaling of the parameters, one can solve the problem in leading order in ϵ.From the solvability condition in the next order one then obtains a closed evolution equation for the cell shape, which due to the cylindrical symmetry can  A1.
The third term on the r.h.s. of Eq. ( 5), P, describes the pushing of actin and can be given explicitly It involves the leading order radial phase field profile ρ0 (ξ) = 1 2 1 − tanh ξ √ 8Dρ and the leading order radial actin polarization field Again, α is a rescaled version of α, hence P ∝ βα is proportional to both the actin nucleation rate, β, and the actin ratcheting velocity, α.
Finally, the last term in Eq. ( 5) describes the effect of ligand-receptor adhesion of the cell on the particle, which is proportional to the adhesion strength κ and where G = G(R, θ) is a function capturing the geometry of the engulfed particle i.e. the effects of the bead phase field Φ.More details can be found in Appendix B and in [39].
While the specific expressions of the contributions look complicated, Eq. ( 5) is a single equation on θ ∈ [0, 2π], which is a substantial reduction compared to two fully 3D equations.Importantly, this equation makes all relevant contributions nicely transparent: The first term is proportional to D ρ , which is the surface tension.Both membrane deformation and the volume constraint impede phagocytosis.In turn, actin pushing -given by integration of the function in Eq. (7) which is peaked at the interface where actin is stimulated by the particle -and adhesion drive phagocytosis.
Using initial and boundary conditions that guarantee a regular and smooth solution, Eq. ( 5) can be solved numerically as shown in Fig. 5B,C.Note that this takes just few seconds using Wolfram Mathematica, compared to several hours running on a GPU for the full 3D model.It should be noted, however, that the asymptotic model can not capture a complete bead internalization and is hence limited to the initial stage of phagocytosis (pedestal phase, cup phase and early wrapping).This is for several technical reasons: first, we assumed a single-valued interface function (cf.Eq. (B.3), where we wrote explicitly r = f (θ, t); an implicit parameterization could be used instead in the future).And second, unlike in the full model, we had to include ϵ in the Gaussians defining the static fields Φ, Ψ, cf.Eqs.(B.1), (B.2), in order to avoid boundary layer problem complications both in time and space.

Conclusion and outlook
We proposed a three-dimensional phase field framework to study the physical part of phagocytosis dynamics.The processes included so far are (receptor-mediated) cellparticle adhesion, inducing actin stimulation by the presence of the bead, which in turn leads to tangential actin polymerization.We have shown that the framework is very versatile and can study geometries/situations that have not been studied yet with other approaches, e.g.phagocytosis of multiple particles and by spreading cells.We also found interesting results, for instance, the study of spheroidal particles revealed that the adhesion dynamics and the actin-ratcheting dynamics have different preferences concerning curvature.Finally, we used a perturbation approach in the thin interface limit to derive a reduced equation for the cell shape R(θ, t), making all important contributions transparent: membrane deformation and volume conservation counteract phagocytosis while actin nucleation and pushing, as well as receptor-ligand adhesion, drive it.
Needless to say that the model still misses most of the biological components and is hence far from being able to describe phagocytosis beyond this simple, physical picture.Nevertheless, as phase field methods have been applied to many different systems already, a large toolbox exists to include relevant processes in the near future, which we like to overview in the following.
Membrane: From the physical side, one should definitely improve the description of the membrane.Due to the high curvatures at the advancing phagocytic tip, bending rigidity constitutes an important contribution counteracting phagocytosis.Even more relevant is to include the fact that at the beginning of phagocyotsis, membrane tension is reduced due to the opening of membrane reservoirs and, much slower, lipid synthesis [40,41].Both explicit surface tension and membrane bending can be included into the phase field approach [42,27], but membrane reservoirs are still hard to model, also in other approaches.
Cytoskeleton: In the approach proposed so far, larger beads stimulate more actin and consequently the uptake is not very sensitive to particle size.While cells can indeed engulf particles of their own size, the actin dynamics is obviously oversimplified.One effect that may be relevant is limited actin resources and its time-delayed transport into the phagocytic cup, amenable to modeling via reaction-diffusion-type additional equations.The closure of phagosomes also needs to account for additional components: since long there is evidence that a contractile activity closes phagosomes [43] in macrophages, which was related to myosin motors; see also Ref. [44] for a recent review.While motors apparently are not crucial always [41], it would be interesting to include them into the framework, which is possible by replacing the actin dynamics, Eq. ( 2) by a more detailed "active gel" model [45].Related to this, the advancing phagocytic front also exerts measurable forces on the bead [46], and elastic effects have been evidenced to constitute a bottleneck for phagocytosis advancement [47].It hence could become relevant to also include elasticity in the phase field approach, as developed recently [18].
Signaling: From the biological side, signaling plays an essential role in phagocytosis [1,5,34].In fact, phagocytosis is highly concerted and regulated, involves the recruitment of many constituents and relies on longer-ranged signaling.For instance, the experimentally observed jump between a receptor-diffusion and an imperfect (either receptor-or actin-related) drift regime of the engulfment dynamics could only be modeled by including a dynamic signaling molecule [28].Signaling processes -modeled by reaction diffusion-like components -can be easily included in the phase-field approach [48,19], similar as done recently in the complementary process of macropinocytosis [13] implementing an activator-inhibitor system.
In conclusion, we hope that the here-presented approach can help paving the way for developing a versatile and more complete physics-based framework to model and shed new light on the complex problem of phagocytosis.Both functions are O(1) in the vicinity of the spherical bead's boundary and decay exponentially otherwise.Note that including ϵ in the exponentials allows to avoid complications with boundary layers both in time and space.This is at the cost of a direct comparison to the full simulations becoming difficult.We now use the standard spherical coordinate system, see Fig. 5A, and assume cylindrical symmetry, i.e., independence of all fields on the azimuthal angle φ.Hence we can write ρ = ρ(r, θ), P(r, θ, t) = pr + q θ.The iso-surface of the cell interface is defined as In order to balance the front dynamics with curvature we impose the following scaling that describes the slow dynamics of a large cell [38], t = ϵ 2 t, f (θ, t) = ϵ −1 R(θ, t), ϵ ≪ 1. (B.4) The transition zone variable is then defined as that will appear in the following analysis.
First, we substitute the scalings of length, time and the parameters, Eq. (B.7) into Eqs.(1), (2).Then we write this system using the transition zone variable, Eq. (B.5), using chain rule when needed.Finally, we apply the asymptotic expansion, Eq. (B.8) and collect terms of the same order of ϵ.
The solutions to the leading order are for the phase field and the actin polarization, respectively.Exact expressions for λ p , and λ q are available in [39].At the next order, O(ϵ), we derive the equation for ρ 1 (ζ).From its solvability condition, we then get a closed evolution equation for the phagocyting cell's interface, R(θ, t), as given already in Eq. ( 5), and discussed in the main manuscript.

Figure 1 :
Figure1: Several phagocytosis scenarios that can be treated within the modeling framework.A1, A2) Early and late stage of a cell (r c = 14) engulfing a small spherical particle (r p = 3).B1, B2) A cell engulfing a particle of size similar to its own (r p = 11 vs. r c = 14).Here, the first phase with the formation of a phagocytic cup is clearly visible, followed by a slower relaxation, where the whole cell body rounds up.Snapshots correspond to t = 50s (A1, B1) and t = 400s (A2, B2).C) A spherocylinder (radius r = 4; cylindrical element of length 4r) is engulfed in a "tip-first" configuration; t = 15s.Sufficiently high actin activity leads to full engulfment after about t = 400s.D) Simultaneous engulfment of two particles (of size r p = 4), positioned at an angle of 90 • with respect to the initial center of the cell; early stage at t = 10s.E) A cell was allowed to spread on a flat substrate (yellow) modeled by an additional static phase field.Then a small sphere was approached.Parameters: A, B, D: α = 4; C: α = 6; others as given in table A1.E: as specified in the main text, see section 4.

Figure 4 :
Figure 4: Simultaneous uptake of two spherical beads by an initially spherical cell.Investigated are three configurations with angles 45 • , 90 • and 180 • with respect to the cell's center.A) shows the normalized contact area per bead relative to the single sphere case (the inset shows the absolute relative contact area per bead).B) shows the normalized overall actin (per particle) as a function of time.45 • is faster, whereas 90 • and 180 • are slower than single bead uptake.C1-C3 show snapshots of the three cases at equal time, t = 50s.Parameters: α = 4, β = 3α, r p = 4, rest as given in tableA1.

CFigure 5 :
Figure 5: A) A schematic description of the geometry specifying the used spherical coordinates, the transition zone variable and the boundary conditions.B) Numerical solution of the reduced model, Eq. (5), parameterizing the phagocyte shape as R(θ, t) for times as indicated (initial condition, very early pedestal phase, maximum engulfment).The spherical particle (radius r 0 = 1) is shown in green and the gridded surface represents the cell of initial radius R 0 = 5.C) Shown are cross-sections corresponding to the stages of engulfment shown in B. Parameters (see also appendix Appendix B): ϵ = 0.1; β = 25, α = 50 (implying α = 5), μ = 8; others as in tableA1.