A model of actin-driven endocytosis explains differences of endocytic motility in budding and fission yeast

A comparative study (Sun et al., 2019) showed that the abundance of proteins at sites of endocytosis in fission and budding yeast is more similar in the two species than previously thought, yet membrane invaginations in fission yeast elongate twofold faster and are nearly twice as long as in budding yeast. Here we use a three-dimensional model of a motile endocytic invagination (Nickaeen et al., 2019) to investigate factors affecting elongation of the invaginations. We found that differences in turgor pressure in the two yeast species can largely explain the paradoxical differences observed experimentally in endocytic motility.


Supplemental Text
Model and results: mathematical details 1. Nucleation and polymerization kinetics and transport of actin.
The model combines detailed kinetics of actin nucleation and polymerization (Berro et al., 2010) and mechanics of actin filament meshwork approximated as a visco-active gel (Kruse et al., 2005). It is formulated in a continuous approximation that does not resolve individual filaments, but rather describes a distribution of actin in the patch by a continuous density of polymerized subunits  as a function of spatial location r and time t , ) , ( t r  . The actin density  is determined by concentrations of all of the species of actin in an actin patch including the newly polymerized ATP-bound subunits (FATP), the subunits aged by ATP hydrolysis and phosphate dissociation (FADP), the subunits bound by cofilin (FCOF), the filaments barbed-ends, both active and capped (BEa and BEc, respectively), and the slowly depolymerizing pointed ends (PE): The concentrations ] [ X are governed by transport-reaction equations of the following type, (S1) The first term in Eq (S1) is the rate of change of ] [ X due to advective flow of actin meshwork characterized by the vector field of local time-dependent velocities ) , ( t r v v  . Like ] [ X and surface densities ] [Y of actin-nucleating species (see below), actin velocities ) , ( t r v are the model variables; they are governed by the mechanics of actin filament meshwork described in the next section. The second term in Eq (S1) is the sum of rates of all reactions affecting a given X . The wiring schematic of reactions involved in nucleation, polymerization, and severing of actin, and capping of the barbed ends of actin filaments is shown in Figure S1. The diagram includes two more volume species that do not directly contribute to ) , ( t r  : the active Arp2/3 complex (ActiveArp) and the Arp2/3 complex in the filamentous actin network (FArp). While the equation governing FArp is of the form of Eq (S1), the equation for ActiveArp has an additional diffusion term: is nonzero only in the vicinity of the rings of nucleation promoting factors (NPFs), as described in the next section (for reasoning and details, see (Nickaeen et al., 2019)).
Reaction steps leading to formation of ActiveArp occur on the surface of the membrane, where they localize to the rings occupied by the nucleation promoting factors WASp ( Figure S1). They involve dimers of WASp bound to G-actin monomers ( WGD ), Arp2/3 ternary complexes consisting of Arp2/3 complex bound to WGD (ArpTernCompl), and activated Arp2/3 ternary complexes (FArpTernCompl), whose surface densities are described by rate equations, where ] [Y is the surface density of a membrane-bound protein Y in molecules per squared micron. Note that while these variables are governed by ordinary differential equations, they also depend on spatial coordinates, given that Y R are nonzero at the rings occupied by WASp, and Reaction rates structure in X R and Y R are as in (Berro et al., 2010). The on (+) and off (-) rate constants of polymerization, capping, cofilin binding, and cofilin-dependent severing from (Berro et al., 2010) were modified by the factor , to reflect the effect of molecular crowding; and the rates of polymerization and capping were modified by an additional factor   ) , to take into account that polymerization and capping slow down under load, where (  is the length of actin subunit; for derivation details, see (Nickaeen et al., 2019)). The expressions of X R and Y R used in computations are shown below: Figure S1. Reaction diagram of nucleation and polymerization of patch actin proposed in (Berro et al., 2010), with added partitioning of species between membrane and cytosol (adapted from (Nickaeen et al., 2019)). Directions of arrows towards or from reaction nodes (yellow squares) determine roles of molecular species (green circles) in a particular reaction as reactants or products, and reactions without products describe disappearance of reactants from the patch. Species connected to reactions by dashed curves act as 'catalysts', i.e. they are not consumed in those reactions.
and the zero subscript denotes constant concentrations and surface density. Table S1 summarizes the values of constants used in the reaction rate expressions and initial conditions for Eqs (1) and for the rate equations describing actin nucleation. For their sources, see Tables 1 and 2 of (Berro et al., 2010). Table S1 also includes the constants involved in actin mechanics: the actin subunit length  and the constants active  and viscous  setting the scales of active and viscous stress and of actin velocities; for derivation of active  and viscous  , see (Nickaeen et al., 2019).

Mechanics of actin meshwork.
We interpret mechanics of actin meshwork as that of a compressible visco-active fluid (Kruse et al., 2005;Prost et al., 2015). In a viscosity-dominated environment of the actin patch, actin velocities ) are governed by the balance of local active and viscous forces (per unit volume), active f and viscous f (Kruse et al., 2005): are the viscous and active stress tensors, respectively ( Î is the unit tensor). In the expression for the viscous stress tensor, is the local average length of actin filaments, computed as  N , with the number of subunits in a filament is the dynamic viscosity of the actin filament meshwork, and ) (  a is the energy (per unit volume) stored in the meshwork during polymerization. By analyzing rheological properties of overlapping actin filaments (Gardel et al., 2003;Kasza et al., Table S1. Model parameters * Maximum surface density of free WASp in the rings is the number of WASp molecules, that corresponds to concentration 0.23 M in a sphere of radius 0.15 m (Berro et al., 2010), divided by area of two rings, 4 , where is the width of one ring. The value is shown for = 5 nm, used in all simulations, and = 3 nm (some simulations were run with neck radii of 6 and 10 nm). ** Upon offsetting by -20 s, becomes -2.8 s (Berro et al., 2010 note that under the filament overlapping condition 1 ) ( 3 / 1 2   N (Doi and Edwards, 1998), the equations (S3) describing actin viscosity crosses from a linear dependence on  over to L 2    (Nickaeen et al., 2019).

Computation domain, boundary and initial conditions.
Equations (S1), (S1*), and (S2) are to be solved in a section of the cell embedding a head-neck invagination, see the diagram of a 3D computational domain  in Figure S2A (drawn not to scale). The domain is bounded by the plasma membrane membrane   , which includes the invagination, and two additional surfaces that define its vertical cylindrical sides and the horizontal floor. In order for the conditions at the outer limits of the domain, which are arbitrary, to not have significant effect on the results, we solved the model in a sufficiently large, 0.5 m in each coordinate direction, neighborhood of the invagination. Note that because the invaginated portion of the domain boundary can move, the shape of the computational domain generally changes with time.
Because the model geometry, localization of membrane-bound species, and corresponding fluxes remain axisymmetric throughout the elongation process, solutions of the model are also axisymmetric. This reduces the problem to solving an equivalent 2D model in cylindrical coordinates ) , ( z r in the domain  depicted Figure S2B (not to scale). The full 3D geometry is then restored for any time by revolving the 2D domain  around the axis of symmetry 0  r .
The equations (S1) for all volume variables, except for [ActiveArp], are of the hyperbolic type and subject to zero-flux boundary conditions at symm_axis on invaginati cell_membr Figure  S2B). For implementation details, see Supplemental Text of (Nickaeen et al., 2019).
Eq (S1*), governing [ActiveArp], is of the parabolic type, because it includes a diffusion term. Therefore, boundary conditions for this equation must be specified at all boundaries in Figure  S2B. Note also the disappearance of the advection term from the boundary condition, which is due to the no-slip boundary condition for Eq (S2) that we describe later in this section.
is defined as follows: where the membrane velocity ) , ( t r u may be nonzero only for the points of on invaginati   . In the model, these points move with the same velocity, were set to zero, and at boundary_2 boundary_1 , the velocities satisfied the zero-stress boundary condition: reflected a small number of seed filaments (Chen and Pollard, 2013).

Computation of the driving force.
The elongation rate is determined in our model by the interplay of the driving force ) ( drive t f produced by the stresses in the assembling actin patch, and the resistive force ) ( resist t f due to turgor pressure (the latter is defined by Eqs (2) and (2*) of the main text). From fluid mechanics, the driving force is found by integrating the viscous and active stress tensors projected on the outward normal (Landau and Lifshitz, 1989), where the integral is carried over the time-dependent invaginated membrane ) ( is the outward unit normal vector to the membrane (directed towards the interior of  ), z e is the unit vector parallel to the axis of symmetry, and ds is the area of an infenitesimal element of the invaginated membrane.

5.
Numerical solution of the model.

Numerical methods
We solved the equivalent 2D model numerically using a moving-mesh solver of COMSOL® Multiphysics (COMSOL Multiphysics, 2015). The solver utilizes Arbitrary Lagrangian Eulerian (ALE) methods based on finite-element (FE) discretization. The ALE methods are described in numerous publications, see e.g. (Donea et al., 2004).
In an ALE simulation, the computational grid points move with the velocities defined as follows. At the moving interfaces, they coincide with the velocities of the points of the interface. In the interior of the domain, the velocities of the grid are arbitrary, so long as they comprise the smooth vector field that would maintain mesh quality throughout the simulation, while preserving mesh connectivity. Correspondingly, the governing equations formulated with respect to a fixed (Eulerian) coordinate system should be reformulated based on the ALE methodology. For the implementation of our model in the ALE framework and details of its FE formulation, see the supplemental material of (Nickaeen et al., 2019).

Nonlinear solver, time stepper, and computational mesh
The coupled nonlinear system of equations describing the cytoplasmic species, Eqs (S1, S*), the membrane-bound species, Eqs (S1**), the ring's position, Eq (1) of the main text, and the actin velocities, Eqs (S2), the mesh update equations, along with the corresponding boundary conditions, were discretized using FE mesh and solved in a segregated manner in COMSOL. Details of such a segregated approach can be found elsewhere (Nickaeen et al., 2017;Ferziger and Perić. 2002). Briefly, the solution of the coupled nonlinear system of equations at each time step is broken into 4 consecutive steps: the Mesh update, the cytoplasmic concentrations update, the membrane ODEs update, and the actin velocity update. At each segregated solver's step, we solve equations for a group of variables while treating the coupling terms consisting of other variables using a fixed-point iteration. Each step of the segregated solver may require additional linearization to account for nonlinearities related to the variables that are solved during that step. In these cases, Newton's method with a constant damping factor of 1 was used and the system's Jacobian was updated at each nonlinear iteration. The resulting linearized system was solved monolithically using a direct MUMPS solver with default solver parameters. Once all steps (one cycle) of the segregated solver is completed, the errors are checked and the iterations for additional, if necessary, cycles are performed to achieve the desired convergence criteria.
Note that even though the force-balance equation does not involve time derivatives, the coupled system constitutes an initial-value problem, so that initial conditions must be specified for all variables (initial values of the actin velocities were set to zero). The time-dependent system was solved using a backward-differentiation time-stepping method of order 1-2. Relative and absolute tolerances of the time-stepper were set to 1 × 10 and 1 × 10 , respectively. Other default solver parameters were used without modification. We verified, by solving the problem with varying solver parameters (including the tolerances of the time-stepper), that the solutions did not depend on specific choices of parameters of the solver. The computational domain was discretized using hybrid quadrilateral and triangular finite elements ( Figure 3A) to accommodate large curvature changes and to allow for sufficient mesh refinement in regions near the invagination in the new head-neck geometry ( Figure 3B).
Full details of the numerical implementation of the model can be found in the attached COMSOL Multiphysics file (Figure2.mph). This simulation file corresponds to the head-neck geometry and is set up with parameters corresponding to Figure 2 of the main text.
6. Effects of the timing of resistance reduction on beginning of elongation and maximum driving force.  Figure S4). As 0 t increases, the point of intersection 'slides' up along the graph of ) ( drive t f , yielding higher begin t and ) ( begin drive t f . However, the effect saturates, once the point of intersection approaches the upper plateau of ) ( resist t f , as the 'sliding' Figure S4. Effects of on and ( ). As increases, the graph of ( ) (dashed curve) shifts to the right (Eq (2) of the main text), and the coordinates of its intersection with the graph of ( ) (solid curve), and , also go up. However, this effect saturates, once the point of intersection reaches the upper plateau of ( ), as is the case with ; is an example of a a 'near-saturation' value of .  Figure S4. From Eq (2)