Role of molecular turnover in dynamic deformation of a three-dimensional cellular membrane

In cells, the molecular constituents of membranes are dynamically turned over by transportation from one membrane to another. This molecular turnover causes the membrane to shrink or expand by sensing the stress state within the cell, changing its morphology. At present, little is known as to how this turnover regulates the dynamic deformation of cellular membranes. In this study, we propose a new physical model by which molecular turnover is coupled with three-dimensional membrane deformation to explore mechanosensing roles of turnover in cellular membrane deformations. In particular, as an example of microscopic machinery, based on a coarse-graining description, we suppose that molecular turnover depends on the local membrane strain. Using the proposed model, we demonstrate computational simulations of a single vesicle. The results show that molecular turnover adaptively facilitates vesicle deformation, owing to its stress dependence; while the vesicle drastically expands in the case with low bending rigidity, it shrinks in that with high bending rigidity. Moreover, localized active tension on the membrane causes cellular migration by driving the directional transport of molecules within the cell. These results illustrate the use of the proposed model as well as the role of turnover in the dynamic deformations of cellular membranes. Electronic supplementary material The online version of this article (doi:10.1007/s10237-017-0920-8) contains supplementary material, which is available to authorized users.


Introduction
Cellular membranes dynamically extend and shrink by removing and replacing their molecular components such as the constituents of their phospholipid and lining cytoskeletons (Staykova et al. 2011). For example, during the lamellipodia formation, the cytomembrane projects outwards as lipid molecules diffuse from within the surrounding membrane and are transported from the cytoplasm (Keren 2011). Moreover, in the process of epithelial mitosis, cells deform to become round as they decrease their surface area, with lipid molecules being transported into the cytoplasm (Raucher and Sheetz 1990). In these cases, membrane deformations are actively driven by forces generated within the cell, i.e., actomyosin contractile forces. Importantly, in general, the membrane morphology is geometrically constrained by its volume-area balance (Ghosh and Singh 1992). Owing to this constraint, deforming cellular membranes requires changing the number of membrane molecules. Thus, molecular turnover has a crucial role in regulating cellular membrane deformations.
This turnover can be regulated by membrane-associated proteins (Peters et al. 2016;McMahon and Gallop 2005;Kozlov et al. 2010). On a molecular scale, these proteins are collectively localized on the membrane and form complexes to pinch off parts of the membrane as liposomes. These liposomes are transported to another membrane and fused through the activities of membrane-associated proteins, such as SNAREs (Grant and Donaldson 2009). Notably, these proteins are known to play the role of mechanosensors, owing to the dependence of adhesion upon the membrane-stress state (Kozlov et al. 2010). Therefore, understanding the effects of turnover upon membrane deformation requires analyzing the feedback from the membrane-stress state on the turnover.
Several computational methods have been proposed to analyze the dynamics of cellular membranes. At a molecular level, molecular dynamics methods have been often used, whereby individual atoms are expressed as particles. Because the membrane dynamics are realistically expressed on the scale of individual lipid molecules, turnover in these methods physically results from molecular interactions (van der Ploeg and Berendsen 1982;Heller et al. 1993;Chiu et al. 1995). On the contrary, at a continuum level, several coarse-graining models have been proposed, ignoring the degrees of freedom of individual lipid molecules (Gompper and Kroll 1997;Ho and Baumgärtner 1990;Boal and Rao 1992;Gompper and Kroll 2004;Zhao and Kindt 2005). In particular, triangular mesh models have often been used to analyze the macroscopic dynamics of organelles and cytomembranes (Noguchi and Gompper 2005b;Ramakrishnan et al. 2013Ramakrishnan et al. , 2015, and also applied to the analyses of cellular mechanotransduction (Atilgan and Sun 2007;Powers et al. 2012Powers et al. , 2014. In these models, membrane morphology is expressed by a triangular meshwork and membrane fluidity is successfully expressed by dynamically remodeling the meshwork topology (Gompper and Kroll 1998). Thus, modeling the membrane turnover on a triangular mesh model yields a powerful tool for investigating the dynamics of cellular membranes.
In this study, we propose a new computational model for simulating the turnover-dependent dynamics of threedimensional cellular membranes. Firstly, we propose topological operations on a triangular network to express the plastic extension and shrinkage of a membrane. Secondly, we propose stochastic descriptions of molecular transport that depend upon the membrane-stress state. Thirdly, using the proposed model, we demonstrate computational simulations of several membrane dynamics and explore the effects of turnover upon membrane deformation. Finally, we discuss the applicability of the proposed model and report new findings on the effects of turnover upon membrane deformations.
2 Multiscale modeling of turnover-dependent membrane dynamics

Description of three-dimensional membrane deformation
The membrane shape is expressed by a triangular meshwork (Fig. 1a), whereby the membrane surface is expressed by a patch of triangles (Fig. 1b, c). The membrane dynamics are governed by an equation for the motion of vertices. By representing the position vector of the ith vertex by r i , the vertex motion obeys the over-dumped Langevin equation as follows: The left-hand side of Eq. (1) indicates a frictional force exerted on the ith vertex. Scalar η is a friction coefficient between the membrane and its microenvironment. The righthand side of Eq. (1) denotes the energetic force acting on the ith vertex, where U is an effective energy function. Variable w i is the Gaussian noise exerted on the ith vertex that satisfies the following statistics: where ... is a statistical average, 0 is the zero vector, and 1 is the second-order unit tensor. Constant k B is the Boltzmann constant and T is the effective temperature. Effective energy U is given by where U cc , U eff , and U act are the constraint, effective, and active energy functions, respectively. From a numerical viewpoint, to maintain a discrete size for the triangular mesh, we describe the constraint energy, U cc , by the following function: where δ [α] is a binary function of the condition α. Here, K r and l rep are the repulsive modulus and distance, respectively. Membrane mechanics are described by the effective energy, U eff . Here, we introduce the ith cell volume v i , the ith triangle area a i , the mean curvature around the ith vertex M i and the mean surface area around the ith vertex A i . Using these variables, U eff is simply described by The first term indicates the volume elastic energy of individual vesicles, where K v and v eqi are the volume elasticity and the equilibrium volume of the ith vesicle, respectively. The second term indicates the surface elastic energy of the membrane exerted on individual triangle areas, where K a and a eq are membrane-surface elasticity and equilibrium area, respectively. The third term indicates the bending rigidity of the membrane, as exerted on individual vertices, where K c is a membrane bending rigidity (Julicher 1996;Tsubota 2014). Variable M i denotes the total mean curvature around the ith vertex: is the jth triangle surrounding the ith vertex.

Description of membrane fluidity and turnover
Cellular membrane has a fluidity, which causes viscous dissipation during membrane deformation. Moreover, membrane molecules are transported between the target vesicle and the reservoir comprising other vesicles within cell. In this model, we regard individual triangular elements as comprising a constant number of membrane molecules. Namely, in the model, the molecules composing the target vesicle is explicitly expressed as the triangular elements, whereas we implicitly express the molecules within the reservoir as a variable. Then, the fluidity can be expressed by a rearrangement of the triangular network (Gompper and Kroll 1998;Gompper 2004, 2005a). Moreover, the turnover can be regarded as a conversion between the triangular elements of the membrane and the number of molecules within the reservoir. Namely, the vesicle size increases when molecules are transferred to and decreases when transferred away from its surface. Therefore, in the model, the fluidity and turnover are expressed from two standpoints: topology and mechanics.

Spatiotemporal scales of our scope
To express membrane fluidity and turnover, we argue spatiotemporal scales of our scope.
From a coarse-graining viewpoint, the target vesicle is composed of the large number of triangular elements. Hence, by representing the current number of triangular elements composing the target vesicle by N t , N t satisfies the following relationship: where N is a positive finite value much smaller than unit. Hence, the spatial scale of vesicle is much larger than that of local triangular elements. Based on the spatial scale, we consider the relationship of the scales of the numbers of molecules in the target vesicle, reservoir, and the triangular element. Here, the constant number of molecules in the individual triangular elements is represented by m u , the total number of molecules composing the target vesicle by M t , and the total number of molecules in the reservoir by M r . Then, M t can be denoted by Because the membrane molecules are transported between the target vesicle and the reservoir, M t and M r are approximately on the same scale: Therefore, M t , M r , and m u satisfy the following relationship: By focusing on the local spatial scale of individual triangular elements, we consider the timescale of membrane fluidity and turnover, represented by τ f and τ t , respectively.
Here, τ f is the characteristic time while the local state of individual triangular elements relaxes by the diffusion of membrane molecules, i.e., the viscous dissipation. τ t is the characteristic time while the local state of individual triangular elements changes by the transportation of membrane molecules from/to the reservoir. Because the transportation of membrane molecules also requires their local diffusion, the time of local fluidity should be much faster than that of local turnover. Otherwise, the membrane should be locally broken during turnover. Therefore, τ f and τ t should satisfy the following relationship: where τ is a positive finit value much smaller than 1. This local relationship corresponds to the global cellular behaviors: While τ f is the period over which the size of the whole cell seems approximately constant, as in blebbing, τ t is the period over which the whole cell size dynamically varies, such as during proliferation or differentiation. Notably, Eq. (11) indicates the timescales of the local dynamics within individual triangular elements, but not the global dynamics of the target vesicle. The timescales of the whole vesicle dynamics caused by the fluidity and turnover, represented by τ F and τ T , can be estimated as respectively. By comparing the amounts of N and τ , from Eqs. (8), (11), (12) and (13), the following relationship is given: In case with τ t /τ f ≥ N t , effects of turnover on membrane are immediately relaxed within the whole vesicle. Hence, in the mesoscopic timescale τ T > t > τ f , the vesicle dynamics seem to be dominated by turnover through the number of membrane molecules. On the other hand, in case with τ t /τ f < N t , effects of local turnover on membrane spend time to relax within the whole vesicle. Hence, in the mesoscopic timescale τ T > t > τ f , the vesicle dynamics seem to be affected by both turnover and fluidity. While the model covers the both cases, this study mainly focuses on the dynamics in case with τ t /τ f < N t , by assuming large vesicles such as cytomembrane.

Expression of membrane fluidity
From a topological viewpoint, the membrane fluidity is modeled by flipping the edges of a couple of neighboring triangles (Gompper and Kroll 1998) (Fig. 2d). Moreover, from Eq. (11), the number of membrane molecules within individual vesicles can be approximately regarded as constant during viscous dissipation. Hence, according to a statistical mechanics, the local state around each edge can be regarded as obeying the canonical ensemble. Thus, the flip frequency of the ith edge, represented by P fi , is given as the following probability: where variable i U indicates a gap in the total energy before and after flipping the ith edge. Notably, τ f reflects the magnitude of membrane viscosity Gompper 2004, 2005a).

Expression of membrane turnover
From a topological viewpoint, the turnover is modeled differently for the increase and decrease in the vesicle size. The increase in the vesicle size is expressed by the local network extension: splitting a couple of neighboring triangles (Fig.  2b). The decrease in the vesicle is expressed by the local network shrinkage: merging a couple of neighboring triangles with their surroundings (Fig. 2c). In these process, the position of a new vertex is determined as the center of the edge shared by the neighboring triangles. Moreover, from Eq. (11), the local state around each edge, from a statistical-mechanical viewpoint, can be regarded as obeying a grand canonical ensemble. Here, we introduce the difference in effective energy before and after the split and merge around the ith edge, represented by i E, as well as the chemical potential of the reservoir, represented by μ r . The frequency of the split and merge around the ith edge, represented by P ti , is given by the following probability: where the sign (∓) is negative for splitting and positive for merging. The factor of 2 multiplying the chemical potential originates from the number of triangles transformed by splitting and merging operations. with surrounding triangles by being altered by the addition of two edges (thick lines in bottom). d Thermodynamic system of local membrane fluidity. e Thermodynamic system of local membrane turnover. In d, e, by regarding the target vesicle as a set of local systems, we consider the systems over a short timescale of fluidity and a long timescale of turnover, respectively 3 Computational simulation of turnover-dependent membrane dynamics

Introducing turnover behavior into the proposed model
In order to analyze effects of turnover on vesicle dynamics, we simply model the mechanosensing regulation using i E and μ r in Eq. (16).
Because the local molecular transport depends on the local stress state of membrane, i E should be a function of vertex positions. As an example, we simply consider the dependency of turnover upon surface-area strain. Because the processes of vesicle fission and fusion require the activities of membrane-associated proteins, i E should involve an active energy cost in addition to passive energy difference such as the change in membrane curvature energy. Moreover, because the density of molecules composing the membrane is likely to be constant, turnover frequency seems approximately proportional to membrane-area strain. Hence, using the first-order approximation, we suppose a linear dependence of turnover upon membrane-area strain. Therefore, by introducing the average surface area of two split or merged triangles adjacent to the ith edge, represented by a i , we define i E as follows: where the sign of (∓) is negative for splitting and positive for merging. Constant t is the energetic cost of molecular turnover. The constant γ t is a critical strain for the energetic reduction of molecular turnover.
To simply define μ r , we consider an isolated system composed of the target vesicle and reservoir. Because membrane transport is actively driven by membrane-associated proteins, the active fluctuation of the turnover is much larger than the thermal fluctuation. Hence, by ignoring osmotic pressure, we focus on the active fluctuation around equilibrium.
For simplification, we assume that the probability density function of M t as the Gaussian distribution with a mean M eq and standard deviation M inst under the condition with i E = 0. Here, M inst can be regarded as the instability of the number of molecules within the target vesicle. This assumption corresponds to defining μ r as follows: where G r is the Gibbs free energy of the reservoir. Here, we introduce the current number of molecules within the vesicle, represented by M t , as a continuum quantity. Then, to convert the triangular elements into molecules, we substitute Eq. (8). By assuming the mass conservation law in the total number of molecules within cell, represented by M tt := M r + M t , M tt is constant. Therefore, μ r can be rewritten as follows: where M req is the equilibrium number of molecules in the reservoir: M req = M tt − M eq . Namely, the employed assumption of M t corresponds to the second-order approximation of the active fluctuation of the number of molecules within the reservoir under equilibrium.

Non-dimensionalization and parameter setting
To solve Eq. (1), parameter values are normalized to have unit length (l), unit time (τ ), unit number of molecules (m), and unit energy (k B T ). Here, l, τ , and m are set as l = a eq 1 2 , τ = 0.1ηa eq /k B T , and m = m u . Hereafter, physical parameters are described as dimensionless values. In case where a specific membrane is focused upon, the physical parameters employed in the simulations can be determined based on those measured by experiments. By assuming the system temperature to be 310 K, unit energy k B T becomes 4.3 × 10 −20 J. Based on this, the values of K c employed in this study correspond to 4.3 × 10 −20 -4.3×10 −19 J. These values have the similar range of those of the dimyristoylphosphatidylcholine (DMPC), the plant thylakoid lipid digalactosyldiacylglycerol (DGDG), and other general lipid membranes (Duwe et al. 1990;Engelhardt et al. 1985;Schneider et al. 1984;Mutz and Helfrich 1990;Evans and Rawicz 1990;Kummrow and Helfrich 1991).
To establish that the proposed model successfully recapitulates turnover-dependent membrane dynamics, several parameters are varied, such as K c and τ t . The state under the initial condition is set as a single vesicle composed of 1000 triangles, which are equilibrated under the condition K c = 10. The equilibrium volume of the vesicle is set to v eq = 2527, which corresponds to the volume of the vesicle with the area 1000a eq and sphericity 0.85. The setting of the geometric constraint l rep in Eq. (5) is described in Appendix A. Moreover, to set physical parameters, the force balance among individual energy terms in Eqs. (5) and (6) is taken into account, as described in Appendix B. Numerical implementation and calculation is described in Appendix C. All model parameters are shown in Table 1.

Proposed model successfully recapitulates turnover-dependent membrane dynamics
To establish whether the proposed model successfully recapitulates turnover-dependent membrane dynamics, we simulate vesicle dynamics in case with and without membrane turnover (τ t = 1.0 and +∞) (Fig. 3). In the case without membrane turnover, the vesicle is slightly deformed by fluctuations while maintaining its surface area (Fig. 3a, c). On the other hand, in the case with turnover, the vesicle is significantly deformed as its surface expands (Fig. 3b, c). Hence, the large deformation is permitted by the surface-area extension (Fig. 3c). Moreover, for a long timescale, the total surface area and number of molecules reached the plateau (Fig. 3c, d). This tendency seems independent on the values of the bending rigidity K c and transport instability M inst ; meanwhile, we could not observe the plateau in the time range of our simulations in case with large M inst . Importantly, the extension is caused by the increase in the number of molecules within the vesicle (Fig. 3d), but not by elastic deformation.
In the process of this expansion, the number of molecules within the vesicle gradually increases in a stochastic manner (Fig. 3e). These results suggested that the proposed model successfully recapitulates the turnover-dependent membrane dynamics.

Molecular turnover adaptively facilitates vesicle deformation
Next, to investigate the effects of molecular turnover upon membrane deformation, we analyze the effects of the bending rigidity K c and the instability modulus of the number of molecules within vesicle M inst (Fig. 4). To focus on the fundamental effects of turnover on membrane dynamics, we set U act = 0 in Eq. (4) by assuming a simple membrane behavior. The vesicle dynamics are found to drastically vary with respect to K c and M inst (Fig. 4a, b). In the case with small M inst , the vesicle slightly deforms while maintaining its surface area. On the other hand, in case with large M inst , the vesicle morphology drastically varies as it extends or shrinks with respect to K c . In case with small K c , the vesicle deforms to be lobate as its surface area increases. In case with large K c , the vesicle deforms to be spherical as its surface area decreases.
To analyze effects of K c on membrane shape, we measured the averaged local Gaussian curvature over every vertex, which is estimated from a set of the surrounding vertices. Notably, the averaged local Gaussian curvature cannot be conserved in the defiance of the Gauss-Bonnet theorem but dynamically vary. This is because the curvature is locally defined at individual vertices, whose number dynamically varies by turnover. Interestingly, the dependence of the averaged local Gaussian curvature on M inst changes directions at three areas; positive in cases K c 2, negative in cases 2 K c 20 and positive 20 K c (Fig. 4b). On the other hand, the number of membrane molecules is inversely proportional to K c independent on M inst . These behaviors can be simply explained by the geometric constraint imposed by the volume-area balance; the surface must be finely folded in cases with large area (K c 2), laminarly flatted in cases with middle area (2 K c 20) and smoothly spherical in cases with small area (20 K c ). Therefore, the resulting membrane morphologies are regulated by the number of membrane molecules through turnover.
The turnover is dependent on K c because it reduces the local residual stress generated by the global force balance: in the case with low K c , the membrane-surface area tends to expand because of the thermal fluctuation force. As K c increases, the membrane-surface area tends to decrease to minimize its bending energy. Therefore, the strain of the membrane-surface area a i /a eq is inversely proportional to the bending rigidity, K c . By sensing local stress as Eq. (17), the molecular turnover is biased to cause expansion or shrinkage. These results suggest that the turnover serves to adaptively facilitate membrane deformation depending upon the membrane-stress state.

Molecular turnover permits autonomous cell migration
Finally, to demonstrate the use of the model, we simulated cellular dynamics driven by an intracellular active force of lining cytoskeleton on membrane. To exert this force, we introduce a locally biased cortical surface energy as follows.
where φ i is the constant defined at each time step as the angle between the x-axis and the vector from the center of the vesicle to the center of the ith triangle. This function is similar to that used in expressing the active energy on cells during collective migration (Sato et al. 2015). As a result, in case without turnover, the velocity of migration drastically decreases (Fig. 5a, c). On the other hand, in case with turnover, the cell dynamically migrates along the x-axis (Fig. 5b, c). This is because the active energy in Eq. (20) drives the transport of molecules from the rear to the front through the turnover (Fig. 5d). This result suggests that cells can migrate even without any traction force on membrane, with turnover playing a crucial role in driving migration. Notably, while the migration can be observed in the wide range of K c and M inst , K c must be high to maintain the spherical cell shape as described in Appendix E.
Notably, because the force from the energy in Eq. (20) is internal within cell but not external, the force is balanced within cell at each time. On the other hand, because φ i dynamically varies with time under membrane deformation, U becomes non-conservative so as to set the system into non-equilibrium. Such non-conservative energy function has been known as to generate active cell movements (Sato et al. 2015). Moreover, the stress-dependency of the turnover in Eq. (16) breaks the detail balance of molecular transport. Thus, despite the force balance within cell, this model can generate cell migration in a physically consistent manner.
In biological systems, there are a lot of types of cell migration such as single and collective cell migration in wound healing, morphogenesis and cancer invasion (Friedl and Wolf 2003). Importantly, while mechanism of the resulting process is a kind of Brownian ratchet, it differs from the mechanism of the well-known single-cell migration. In this mechanism, the front extension of the cell initiates migration: It extends ahead by cytoskeletal polymerization and the rear follows it by cytoskeletal contraction. On the contrary, in our model, localized active tension at the rear initiates migration by directionally transporting molecules to the front. Notably, our model does not conflict with the mechanism of the single- cell migration. Therefore, turnover may contribute to cellular migration as a main or secondary process.

Discussion
In this study, by making several physical assumptions, we integrally formulated the fluidity and turnover of a cellular membrane in a physically consistent manner. In the computational simulations, the cellular vesicles actively deform by the turnover of the membrane molecules (Fig. 3). During deformations, the triangular property was maintained to satisfy the modeling assumptions as in Appendix D. Moreover, based on the dependence of turnover upon the stress state, vesicle morphology drastically differed with the bending rigidity (Fig. 4); a smooth sphere was obtained under high rigidity, and a lobate vesicle was obtained under low rigidity. Furthermore, based on directional molecular transport by turnover, the localized active tension on the membrane drove the cellular migration of the vesicle (Fig. 5). From these results, the proposed model successfully recapitulated the turnover-dependent dynamics of three-dimensional cellular membranes.
In the computational simulations, the turnover behavior described in Eq. (16) was simply provided under the following assumptions: (i) that the energetic cost for molecular transport, i E, depended on the local strain of the membrane area, as shown in Eq. (17) and (ii) that the chemical potential μ t stabilized the number of molecules in the reservoir, as shown in Eq. (18). In general, these assumptions were not necessary for the proposed model. We emphasize that the functions of i E and μ r can be entirely arbitrary.
Although the cellular mechanical property was expressed simply in these simulations, as described in Eq. (6), it can be expressed in more detail. For example, in cells with rich lining cytoskeletons, cell membranes could have non-negligible longitudinal and transverse elasticity, which may be important for their deformations (Heinrich et al. 2001). Moreover, defined as that between the x-axis and the vector from the center of the vesicle to the event site. In c, the frequency density is calculated as the local frequency per the surface area of sphere within the scope angle. In d, the local surface strain is calculated as the strain of individual triangular elements, and the solid line and bar indicate the average and standard deviation within the scope angle, respectively. These dynamics are calculated under the condition K c = 10, M inst = 10, and K act = 3.0 during blebbing, the actin cytoskeleton lining membrane is locally broken, which causes a spatiotemporal inhomogeneity in membrane rigidity (Manakova et al. 2016). These details can be reflected by choosing a proper effective energy function, U eff in Eq. (6). In principle, detailed expressions for turnover at a molecular scale are limited in the proposed model; for example, the adhesion process of membrane-associated proteins cannot be directly expressed. Such behaviors at the molecular scale must be instead coarse-grained into those at the minimum scale of the proposed model. For example, the dependence of protein adhesion upon lipid behavior can be expressed by i E as a dependence of the turnover upon membrane strain. Moreover, the biased adhesion of proteins to lipids can be expressed by μ r as directional molecular transport. Thus, by designing i E and μ r , the proposed model can be applied to various behaviors of turnover.
The proposed model will help researchers understand the mechanics of 3D cellular dynamics involving membrane turnover, particularly the causal relationships between single-cell dynamics and the underlying molecular transport. The proposed model is a powerful approach for addressing the manner in which molecular turnover affects 3D cellular dynamics at a subcellular scale. Understanding these functions of turnover is necessary for better understanding cellular dynamics, and this will be useful as fundamental knowledge for controlling these dynamics in bioengineering. Moreover, by combining the proposed model at a subcellular scale with the models at a multicellular scale (Okuda et al. 2013(Okuda et al. , 2015, it will be possible to predict comprehensive cellular dynamics from molecules to tissues. Therefore, the proposed model will contribute to exploring the exploration of the frontiers of cellular biomechanics.

Conclusion
To analyze the effect of molecular turnover on threedimensional deformations of cellular membranes, we proposed a new computational model for simulating the turnoverdependent dynamics of three-dimensional cellular membranes. The proposed model successfully simulate turnoverdependent membrane dynamics, and suggested the roles of turnover to drive the adaptive deformation and directional migration of vesicles. These results illustrate the importance of turnover in the dynamic deformations of cellular membranes in addition to the use of the proposed model for exploring general effects of molecular turnover on cellular dynamics.  All experiments were performed on a cluster computer comprising 12 nodes with 2.9 GHz Intel Xeon dual processors and 64 GB RAM.

D Validation of assumptions
In the proposed model, we made several assumptions: (i) that the timescales of membrane fluidity and turnover are related as in Eq. (11); (ii) that there is topological flexibility for membrane fluidity, and (iii) that the force balance is given by Eq. (22). Preassumption (i) is satisfied because it is given by the set of parameters. On the other hand, preassumptions (ii) and (iii) are non-trivial because they result from the dynamic process of membrane deformation. To ascertain whether these preassumptions are satisfied in the results, we analyze properties of the membrane and triangles. Firstly, to test the membrane fluidity, we analyze the probability distribution of the number of sides, which is broad from 4 to 9 (Fig. 7a). Hence, in the case where l rep is set as Eq. (21), the constraint energy permits the topological flexibility of the meshwork. Secondly, to test the balance of the individual energetic forces, we analyze the probability distributions of vesicle volume and triangular area. These distributions approximately equal those determined by the first and second terms of Eq. (6) (Fig. 7b, c), respectively. Hence, the preassumed force balance of Eq. (22) is satisfied. Moreover, from the distribution of triangular area in Fig. 7c, the discrete size of the triangular meshwork is approximately maintained. These results warrant certain predictions of the model from a physical viewpoint.

E Effect of membrane turnover on cell morphology during migration
To analyze the effect of molecular turnover on migration, we simulated cellular dynamics driven by the intracellular active force, described by Eq. (20). Figure 8 shows the cell morphologies during migration with respect to K c and M inst . The results suggested that K c must be high to maintain the cell shape to be spherical. Otherwise, the cell surface strongly undulated.