Cell Invasion Dynamics into a Three Dimensional Extracellular Matrix Fibre Network

The dynamics of filopodia interacting with the surrounding extracellular matrix (ECM) play a key role in various cell-ECM interactions, but their mechanisms of interaction with the ECM in 3D environment remain poorly understood. Based on first principles, here we construct an individual-based, force-based computational model integrating four modules of 1) filopodia penetration dynamics; 2) intracellular mechanics of cellular and nuclear membranes, contractile actin stress fibers, and focal adhesion dynamics; 3) structural mechanics of ECM fiber networks; and 4) reaction-diffusion mass transfers of seven biochemical concentrations in related with chemotaxis, proteolysis, haptotaxis, and degradation in ECM to predict dynamic behaviors of filopodia that penetrate into a 3D ECM fiber network. The tip of each filopodium crawls along ECM fibers, tugs the surrounding fibers, and contracts or retracts depending on the strength of the binding and the ECM stiffness and pore size. This filopodium-ECM interaction is modeled as a stochastic process based on binding kinetics between integrins along the filopodial shaft and the ligands on the surrounding ECM fibers. This filopodia stochastic model is integrated into migratory dynamics of a whole cell in order to predict the cell invasion into 3D ECM in response to chemotaxis, haptotaxis, and durotaxis cues. Predicted average filopodia speed and that of the cell membrane advance agreed with experiments of 3D HUVEC migration at r2 > 0.95 for diverse ECMs with different pore sizes and stiffness.


Introduction
Cell migration in the three dimensional extracellular matrix (ECM) plays a crucial role in a wide variety of biophysical processes, such as wound healing, morphogenesis, angiogenesis, tumor growth, and cancer invasion [1,2]. Migration dynamics in 3D are significantly different from those observed on 2D ECM surfaces [3]. Cells invade into the ECM, extend filopodia into the gel, and degrade and remodel the surrounding ECM. Cells sense the direction and magnitude of complex cues and exhibit multifaceted responses, such as chemotaxis, haptotaxis and durotaxis responses from the 3-D extracellular environment.
In contrast to their behavior on 2D surfaces, filopodia in 3D must penetrate through and interact with the surrounding ECM fibers, probing for an open space or gap in the ECM fibers to extend. When they encounter an ECM fiber, they bind to it via FCs, and subsequently generate traction forces to pull the cell membrane forward. This complex process and the intricate interactions that occur between the filopod and the ECM fibers remain poorly understood. How are these multifaceted activities coordinated to allow first the filopod and eventually the whole cell to penetrate into the 3D ECM? What are the essential features of 2D filopodia dynamics needed to predict 3D migration? To address these questions, we have conducted experiments to probe the detailed interactions as a single filopodium penetrates into a 3D gel matrix, and constructed a computational model to predict coordinated filopodia behaviors and cell-ECM interactions in 3D.
Our experimental observations using 2D and 3D time-lapse data revealed that filopodia a) crawl or slide along ECM fibers near the tip, b) tug on the fibers to generate local forces, causing deformations within the surrounding gel, and c) probe the local fiber network and coordinate their multiple activities: protrusive outgrowth, retraction, and contraction. To explain these behaviors we have built a computational model in which the crawling/sliding of filopodial tips along fibers is described as a continuous process of forming and rupturing focal complexes (FCs) between integrins distributed along the filopodial shaft and ligands on the ECM fibers. The strength of FCs is determined by binding kinetics and other factors, including the local stiffness and porosity of the ECM. Each filopodium alternates between different modes of behavior depending on the properties of FCs. We have established simple rules to describe the switching and coordination of multifaceted filopodia behaviors. The model reproduced not only our experimental data, but also many of the previously reported behaviors, including frictional slippage, and load-and-fail phenomena at hard and soft ECMs [4].
Predicting the filopodia-ECM interactions and the resultant penetration of the whole cell into a 3D gel matrix requires a detailed model of the ECM. We have constructed a network model of cross-linked ECM fibers, with specified pore size and fiber stiffness to match the experimental measurements of bulk elastic moduli. This fiber network consists of straight fibers having both extensional and bending stiffnesses [12,13]. The ECM fiber network degrades due to matrix metalloproteinases (MMPs) [14,15]; the expression of proteolytic enzymes at the cell membrane [16] dissolves the fiber, ruptures the fiber, or ruptures the crosslink [17,18]. Over the ECM network, chemo-attractants diffuse and react to the cell. Furthermore, the mechanics of the cell membrane and nucleus as well as acto-myosin contraction and actin stress fiber formation are all coupled dynamics inside the cell and the filopodium that must be integrated with the filopodial / ECM dynamics.
Previously, many modeling approaches in the areas of cell migration have already provided insights into processes of cell-adhesion and cancer cell invasion. Such approaches include a coarse-grained Langevin dynamics model of lamellipodium protrusion by actin polymerization [19], an invadopodia penetration model with the effect of crosslink on ECM degradation [17], a force-based, individual-based modeling framework that links single cell migration and ECM fibers through contract guidance and matrix remodeling [18], a multi-scale model of dynamic of cell colonies [20], and a cell invasion model in fiber networks and confined microchannels using extended cellular Potts model (CPM) [21].
Our modeling approach exhibits both mechanistic and chemical interactions of 3-D ECM with filopodial and cellular membrane structures since the degradation of ECM is more essential as the cell penetrates into stiffer and denser ECM [22]. Mechanistic interactions of ECM induce local compaction [23], migration and remodeling of individual ECM fibers [24]. On the other hand, chemical interactions of ECM guide dynamic variations of filopodial orientations due to growth factor gradients in the ECM gel as well as the degradation of ECM fibers due to the secretion of focalized proteolysis at the cellular membrane [25]. As further evidence for the model, we also find the densification of ECM fibers in surroundings of the filopodia at both simulation and experiment since filopodia create a substantial traction force and corresponding ECM fibers are considerably deformed. In addition, our modeling approach is developments of previous our modeling works of intracellular mechanics; an integrative cell migration model incorporating FA dynamics, cytoskeleton and nucleus remodeling, actin motor activity, and lamellipodia protrusion was developed for predicting 1) cell migration behaviors on 3D curved surfaces, such as cylindrical lumens in the 3D ECM [26], and 2) cell spreading and migration behaviors on micropatterns and planar substrates with various fibronectin coating concentrations [27].
Integrating these we predicted the average speeds at which the filopodium and cell penetrate into the ECM, obtaining excellent agreement (r 2 = 0.995) with experiment data for diverse pore sizes and gel moduli. To our knowledge, no cell invasive model that takes into account penetration of the cell into a 3D ECM environment has previously been reported that reflects both the cell-ECM and filopodia-ECM interactions.

Cell invasion dynamics model
To create the overall computational model, we integrate four modules, each capturing a different physical aspect influencing migration: 1) filopodia penetration dynamics [4] (Fig 1A); 2) intracellular mechanics, including formation of FAs and actin stress fibers (SFs), and remodelling of cellular and nuclear membranes [26,27] (Fig 1B); 3) reaction-diffusion mass transfer [14,15] (Fig 1C); and 4) dynamics of ECM fiber networks [12,13] (Fig 1D). In particular, it should be noted that FAs are different from FCs in that FAs (1-5 μm in size) are formed on the cellular membrane with a long turnover time > 5 minutes, but FCs (~0.5 μm in size) are  E). B) cellular membrane mechanics showing the membrane is not only connected by actin stress fibers, but also anchored to elastic ECM fibers by forming focal adhesions (FAs) (See 'b' in e), and viscoelastic behaviors in cellular membrane is modeled using Kelvin-Voigt model. C) Schematic diagram of MMP-2 activation [14] and selected simulation results of VEGF, MMP-2 and TIMP-2 concentration distributions using simplified model for MMP-2 activation. D) a magnified schematic of 'c' in E) showing the elastic ECM fiber networks is organized with collagen fibers, crosslinking molecules and free (or noncrosslinked) molecules. E) An integrated schematic representation of MT1-MMP and TIMP-2 secretions at the membrane near the roots of filopodia; directions of outgrowing filopodia are guided when tips of filopodia sense local gradients of VEGF toward a chemotactic cue. formed on the filopodial membrane with a short turnover time < 5 minutes [28]. To incorporate viscoelastic behaviors in cellular membrane, line elements of actin cortex in the cellular membrane can be modeled using Kelvin-Voigt model (a spring and a dashpot together in parallel) (Fig 1B). The detailed equations that govern each of these dynamical processes are summarised in Table 1, and the list of simulation parameters are also summarised in Table 2.

Filopodia penetration dynamics
We identified and modelled four phases of filopodial dynamics. Three resemble the phenomena previously reported: an outgrowing phase due to protrusive polymerization of actin [5][6][7]; a retractile phase due to zero or weak FC forces at the filopodial tip and fast myosin motor activities along the filopodial shaft [9]; and a contractile phase due to strong FC forces at the filopodial tip and slow myosin motor activities along the filopodial shaft [4,9] (Fig 1A). The fourth phase, which we refer to as the "tugging phase", begins with the formation of FCs near the tip of the filopodium. The point of attachment between a filopodial tip and a nearby ECM fiber migrates along the fiber, a tension builds up between them as the filopodium moves, and the bond either ruptures or results in the generation of a significant traction force, depending The four modules describe the dynamics of the filopodia (F), cell (C) and ECM (E) modules and the Reaction-diffusion (RD) modules. F module: i and F f AM;i are the elastic force, focal complex force, force due to protrusive polymerization of actin filaments; and contractile force due to actomyosin (AM) motor activity at the i-th filopodial node, respectively. C module: C is composed of three sub-modules representing the cellular membrane (CC), the transduce layer (CT) and the nuclear membrane dynamics (CN). CC module:F c FA;i , F c E;i , F c L;i , and F c T;i are the focal adhesion force, the elastic energy force, lamellipodium force, and cortical tension force at the i-th cell membrane node, respectively. CT:F t E;i ,F t SF;i , and F t T;i are the elastic force, SF force, and cortical tension force at the i-th force transduce node, respectively. CN: F n E;i and F n SF;i are the elastic force and SF force at the i-th nuclear membrane node, respectively. E: F e FA;i ,F e FC;i and F e E;i are the FA force, FC force and elastic force at the i-th ECM fiber node, respectively. Here, C f , C c , C t , C n and C e are friction coefficients associated with the energy dissipation at F, CC, CT, CN and E modules, respectively. In addition, C cort is a drag coefficient associated with viscoelastic behaviours in the actin cortex. All dynamic models of simulations are carried out using a fourth order Rosenbrock method based on an adaptive time-stepping technique for integrating ordinary differential equations with convergence criterion < 10 −4 . F module is geometrically coupled with CC module and with E module through, CT module is coupled with both CN modules with an equation, F t SF;i þ F n SF;j ¼ 0. RD module is divided into two subgroups according to whether or not the diffusion terms are included (See also Table 2). RD with diffusion: X 1 is the concentration of VEGF, MMP-2 and TIMP-2, or ligand at both the ECM and medium domains. Here, medium domain means the microfluidic channel or blood vessel, and the cell adhere to the interface between ECM and medium domains (See Fig 1C).  on the strength and spatiotemporal properties of the FC formation. This phase plays a critical role in switching among the other phases and coordination of the diverse filopodial dynamics, leading to either success or failure of cell migration depending on local ECM conditions. First, in the outgrowing phase, polymerization of actin filaments generates a force F f P;i on the filopodial membrane [29]. The average magnitude of the maximum protrusive actin pressure is in the range of 5-10 nN/μm 2 . It has been reported that a few nN of force can exerted by a few hundred actin filaments per μm 2 , which implies that each filament contributes on the order of 10-20 pN [30]. Thereby, the imposed maximum magnitude of F f P;i~2 nN (F f P;max ) is reasonable since the diameter of single filopodium is assumed to be 300 nm consisting of >30 actin filaments. It should be noted that F f P;i is exerted only at the tip of the filopodium. To incorporate a chemotactic response from the 3D ECM environment, the direction of F f P;i is predicted from the gradient of the VEGF concentration (C VEGF ) at the tip of the filopodium. Thus,  (Fig 1E).
Second, the contractile phase arises from the acto-myosin contraction. At contractile phase, we assume that a bundle of actin microfilaments is assembled by actin-myosin II interactions [31]. It has been observed that veil (or lamellae) advances following filopodial pulling (or contractile) motion, a process in which two kinds of adhesion mechanisms with collagen fibers in the ECM take place, i.e., FCs formation at the tip of filopodial shaft and FAs formation at the root of filopodium [32]. These two kinds of adhesions at both tip and root of the filopodium are required for gaining stable traction forces from the surrounding collagen fibers in the ECM, and for transmission of the force to allow the advancement of the cellular membrane via the filopodial contractile motion. Here, stable traction forces indicate sufficiently large magnitude to withstand the acto-myosin contractile force. In addition, for the condition of the veil advancement, the focal complex force at the filopodial tip must be stronger than the focal adhesion force at the veil. Otherwise, weak complex force leads to the retractile phase. Therefore, we model filopodia penetration dynamics at the contractile phase by taking into account the complex myosin motor activity and FCs formation at the tip of filopodial shaft.
Third, the retractile phase is also arises from the acto-myosin contraction in a similar manner of the contractile phase. However, the retractile phase is induced by weak traction forces from the surrounding ECM, which result in fast, oscillatory 'load-and-fail' traction dynamics during the retractile phase (S1 Fig and S1 Video). That is, the reconnection of FCs at the tip allows a filopodium to repeatedly probe the surrounding ECM fibers, but the filopodium retracts immediately if it experiences weak traction forces at the tip [10].
Finally, we have found through experiments that there exists another phase in filopodial dynamics: the tugging phase (S2 Video). The tip of a filopodium apparently crawls or slides along a nearby ECM fiber or multiple fibers (Fig 2A). We have observed that the binding site of FCs moves along the ECM fibers and that the bound ECM fibers are pushed or pulled by the filopodial tip. The adhesive force of FCs is sufficient to prevent the filopodium from retracting. This crawling or sliding can be viewed as a continuous process during which FCs form and rupture as they move along the ECM fibers, as depicted in Fig 2B.

Filopodia focal complex formation
Here, we build a computational model to predict the continuous formation and rupture of FCs along the ECM fibers based on binding kinetics between integrins on the filopodial membrane and collagen molecules in ECM fibers. This technique is similar to the one used to predict the formation of a FA at a cellular membrane that interacts with ligands in the ECM fibers [27]. The FC force at the k-th filopodial node is given by where n k b is the number of integrin-collagen bonds, κ LR is the spring constant of a single ligand-receptor bond (~1 pN/nm) [33], L b is the average stretched length of the ligand-receptor bonds, λ is an unstressed length of bonds (~30nm) [34] andn f R;k is a unit vector at the local surface of the k-th filopodial node toward the bonding site between the j-th and j+1-th nodes in the i-th fiber ( Fig 2B). (L b -λ) represents the stretched distance from the equilibrium. From Fig  2B, this intersection position, that is, the root location of receptor and ligand bonds (x L,i ) between x e ij and x e ijþ1 , is given by where h p is the gap between the k-th filopodial membrane node and ECM fiber node, andn w is a vector normal to the curved surface of the cylindrical fiber, defined by ðx e ijþ1 À x e ij Þ Â ððx f k À x e ij Þ Â ðx e ijþ1 À x e ij ÞÞ and shown in Fig 2B. We utilize Bell's model to simulate the stochastic nature of bond rupture and formation. Bell's equation for the kinetic dissociation rate is off is the kinetic dissociation rate (1 s -1 ) under unstressed condition with an equilibrium distance λ, x b is the separation distance between the bound state and the transition state (0.02 nm), k b is the Boltzmann constant, and T is the absolute temperature [35]. In addition, represents the maximum rupturing force exerted on single molecule of ligand-receptor bond (~200 pN). During the dynamic process, the filopodium switches between these four phases. Depending on the strength of the FC force, the length of the filopodium, the duration of the current phase, and other factors, transitions between phases are determined. This can be modelled as a discrete state transition network detailed in S2 Fig. Briefly, the outgrowing phase switches to the tugging phase when FCs are formed within a specified time, otherwise it switches to the retraction phase. In addition, the formation of FCs is assumed to occur when h p is less than 100 nm, and is described by a stochastic process due to binding kinetics between receptors and ligands on the surface of ECM fiber. Monte Carlo simulation methods have been established for various ligandreceptor binding kinetics in the literature [27]. Each focal FC consists of a bundle of ligand-receptor bonds (Fig 2B), each of which ruptures and binds stochastically. Let P b be the probability with which a single receptor binds to a ligand on the ECM fiber during a time interval Δt.
where k f is the forward reaction rate (1 molecule −1 s −1 ), C b represents the density of bound ligands, C L the original density of the ligands (molecules area −1 ), and A L is the local area of a single fiber associated with the integrin node under consideration. Note that (C L -C b ) represents the number of unbound ligands available for bonding in the vicinity of the integrin node. In simulations, values of A L C L were identically set to be 300 molecules per a ECM fiber node in three ECM fibers network models. Once in the tugging phase, the strength of the FCs is tested (rupture test), and the phase switches to the contractile phase if the bond does not rupture and the tension rises beyond a specified threshold level at the filopodial tip (3-4 nN/μm) [36]. If it ruptures, it switches back to the outgrowing phase. The formation of FCs is restricted to the proximal tip of the filopodia, since it is known that veil advance typically results from the formation of FCs proximal to the tip of the filopodial shaft [32] (Fig 1A and 1E, S3 Video).

Geometrical model of contractile filopodia
The filopodial model is geometrically composed of N AM compartments of acto-myosin (AM) assemblies; the first compartment is attached to the root of filopodium and the last compartment is connected to the tip of filopodium (Fig 3A and 3B). We model filopodial contractile motion in a manner characterized by that of actin stress fibers (SF) [27]. The stiffness of an AM assembly is variable, AMs [37], A AM is the average cross-sectional area of AMs in a filopodium and L 1 AM;j is an unstressed length of a single compartment of the j-th AM. The length of each compartment contracts at both ends according to the myosin II sliding rate, v m,j . Therefore, Furthermore, it is known that myosin motors slides on actin filaments in opposite directions of FC forces at the filopodial tip [4]. That is, an increasing elastic load from the ECM is transmitted through FCs at the tip of the filopodium, through the AM assembly, to the myosin motors, which, in turn, decreases the myosin speed (v m ). It has been known that myosin sliding velocities affected by ECM stiffness, that is, myosin sliding speed is high on the soft gel but it is low on the hard gel [4]. In addition, measurements on few or many myosin molecules by optical trap [38] reported that the force-velocity relationship for these molecules is significantly similar to that of whole muscle [39] or muscle fibers [40,41]. To incorporate these characteristics into the filopodia dynamics, we adopt the force-velocity relationship for the whole muscle [39] into the following equation: v m ¼ v m0 F stall ÀF TR F stall þc m F TR where v m0 is the sliding rate of myosin in the absence of load (10 nm/s) [42], F stall is the stall force of 1 nN, c m is a dimensionless myosin parameter of 0.1, and F TR is the sensed elastic force from the ECM at the tip of filopodium. It should be noted that the speed of myosin become zero when F TR exceeds F stall , and linear relationship of force-velocity can be tuned to that of the nonlinear relationship by increasing the value of c m more than 1.
The total elastic energy stored in the AMs in the filopodium is given by where d AM,j represents the distance of the j-th contractile AM compartment under tension ( Fig  3B). Using the virtual work theory, forces due to contractile myosin motor activity at the j-th node of filopodial shaft is given by Characterization of three ECM fiber network models. It is known that mechanical properties of collagen gels depend on the pH of the collagen solution during polymerization [43]: gels become stiffer and structures of gels were found to contain both decreased fiber diameters and pore sizes as the polymerization pH is increased. Accordingly, we built three different models of ECM fiber networks comprising of three different pore sizes of 0.5, 1.0 and 1.5 μm and three different fiber diameters of 28, 34 and 41 nm, respectively (S3 Fig and S4 Video). First, we aimed to simulate mechanical stretching tests for the three ECM fiber network models in order to characterize and compare the three models with experimentally measured bulk moduli of ECM gels [43]. Various simulations of mechanical stretching tests for each ECM fiber network model were performed using two parameters of single fiber diameters and moduli (S4-S6 Figs and S1 Text). We found that the moduli of the model are linearly and significantly increased as single fiber diameter or single fiber modulus is increased in each ECM model (S3A and S3B  Fig). Both simulation and experiment show an excellent agreement over both ECM fiber diameters of 28, 34 and 41 nm and ECM pore sizes of 0.5, 1.0 and 1.5 μm. Statistical analysis of linear regression was performed by comparing the experiment and the simulation in terms of the mean values of modulus for the same conditions of fiber diameter and ECM pore size. Good correlations were found between the two with R 2 = 0.898 (S3C Fig). It should be noted that identical fiber modulus of 1 MPa for all the ECM models gives the best agreement with the experiments. Moduli of ECM gels are strongly dependent on single fiber diameter and ECM pore size, and the above fiber network model successfully verifies these properties.
In literatures, measured single fiber moduli and diameters in aqueous media have been reported to be from 32 to 800 MPa, and from 40 to 1000 nm, respectively [44][45][46]. Interestingly, Graham et al. showed a peak value of 32 MPa at the larger strains than 4 and the lowest strains~2 MPa at the range of strains between 0 and 2 using a fiber diameter of 40 nm [44], and our simulated fiber diameters and modulus are in a good agreement with their experimental data at the low strains <1 since we assume the cell can generate traction force interacting with fibers at this low strain regime. In addition, an image-based multiscale modeling technique used fitting parameters of fiber modulus of 6.7 MPa and a fiber diameter of 100 nm to predict both tissue-level and network-level fiber reorganization [47].
Prediction of 3D cell-ECM interactions and penetration speeds and comparison with experimental observations. Following successful characterizations of the mechanical properties of the three ECM models with collagen gel experiments [43], simulations of filopodia penetration and cell invasion were performed using these three ECM models with three different pore sizes of 0.5, 1.0 and 1.5 μm and three different fiber diameters of 28, 34 and 41 nm ( Fig  4A, 4B and 4C, S5-S7 Videos) to predict coordinated filopodia behaviours and cell-ECM interactions in 3D. This allowed us to investigate the dynamic interplay between filopodial traction generation and ECM gel remodelling in 3D in a quantitative manner.
Furthermore, for comparison and validation of the model, experiments were performed using HUVECs transfected with cytosolic GFP migrating in the three types of collagen gels with pH levels of 9, 7, and 5 ( Fig 4D, 4E and 4F, S8-S10 Videos) filled in-vitro microfluidic device [48]. Penetration speeds were measured at both tip and root of each filopodium located on the leading edge of the cell. Each experiment was recorded at a time-interval of five minutes over 30-60 minutes, while simulations were conducted with a time-interval of one second over 10-20 minutes for the three gel types (Fig 4A, 4B and 4C). For an example, the time-averaged speed at the filopodial tip,v tip mean , is expressed as following: where t 0 is the initial time, t f is the final time, and N sample is the number of samples.x tip i ,y tip i , and z tip i are coordinates of filopodium at the i-th time-interval. We have found that speeds of both tip and root of filopodia increase as the pore size is increased in experiments (Fig 4H). The simulated speeds of both tip and root of filopodia, too, shows a trend similar to the experiments (Fig 4G). The simulations and experiments have shown an excellent agreement for all the ECM pore sizes of 0.5, 1.0 and 1.5 μm (or pH levels of 9, 7, and 5). Statistical analysis of linear regression was performed for comparing the experiments and simulations in terms of the mean filopodial tip speed as well as of the mean root speed under the same conditions of ECM pore size (or pH level). High correlations were found between the two, with R 2 = 0.995 and a slope of 0.7 (Fig 4I).
Filopodia state dynamically alters during penetration. Experimental observation revealed additional behaviours of filopodia that support the mechanisms identified in our filopodia dynamic model as well as those of prior works [4][5][6][7][8][9][10]. Filopodial protrusive and contractile motions were recorded simultaneously with ECM fibers deformation and remodelling ( Fig  5A and S11 Video). An example can be seen in Fig 5A: two filopodia extending in different directions and their proximal ECM fibers were tracked over time: 'Filos A' and 'Filos B'. Filo A started crawling on collagen fibers at time 0, its tip was further protruded and branched along multiple collagen fibers in different directions (2 min, tugging phase); collagen fibers were stretched due to the contractile motion of filopodia (4 min, contractile phase); following that, fluctuating motions at the filopodial tip was further observed (6 min-10 min, tugging phase).
Movements of the filopodium and those of a few specific ECM fibers near the filopodium were highly correlated, indicating the existence of linkage between the two. To investigate this correlation in movements, the background fiber images were eliminated from the original 3D confocal images and displacements of both filopodium tip and neighbouring fibers were tracked over time (Fig 5B and S12 Video). Fig 5C shows how the speed of this ECM fiber marked with a blue sphere is completely synchronized with the speed of the filopodial tip for the first 30 min. Note that plus and minus signs represent forward and backward movements of filopodium. In fact the filopodial tip crawled along this particular fiber; tugging the fiber bit by bit (t+2, t+6, and t+8 min) induced the fiber motion that was completely out of phase with the filopodium (Fig 5C). Interestingly, at t+10 min, we can notice relaxation of the fiber (blue sphere displaces away from the filopodium tip), which is correlated with the filopodium tip switching from that fiber to an adjacent fiber (t+10 min, bottom yellow arrow).
Filopodia state changes were also observed in these experimental data. The filopodium was rapidly retracted at t+18 min (S12 Video and Fig 5D), and fast oscillatory 'load-and-fail' traction dynamics were observed at t+32 min (blue arrows in Fig 5E). Plotting the tip speed together with the speed of the filopodium root, which is approximately the same as the speed of veil (cell membrane advance) ( Fig 5E); we can find how filopodia dynamics influences the cell advance. In addition, Filo B also showed similar behaviours except for the rapid retraction of Filo A (Fig 5F). The migratory trajectory of the cell indicated that Filo B was at the leading edge of the cell and had been polarised in that direction. As a result, its life span appears to be longer, exhibiting slower contractile motion rather than rapid retraction as in the case of Filo A (Fig 5E and 5F).
Comparisons of traction fields in ECM fiber network. Distinct behaviours of filopodial motion have been confirmed using a 3D cell migration microfluidic assay in which 200 nmsized fluorescent beads are embedded into the collagen gel. As a result, beads move towards the extending filopodium tip during the tugging and contractile phases, and away from the filopodium during retractile motion (S13 and S14 Videos). Based on these experimental observations, computational model of filopodia penetration dynamics is validated by showing similar motions of ECM fibers toward or away from the filopodium during the tugging and contractile phases or retractile phase, respectively. These imply that ECM fibers repeatedly underwent tension and relaxation; the direction of displacement vectors changed accordingly, switching between pointing towards the filopodia and pointing away from it cyclically (Fig 6A and S12  Fig and S15 Video). Furthermore, the number of receptor-ligand bonds at the tip of filopodium, the length of filopodium, and the strength of traction force at the tip of filopodium were computationally measured, as shown in Fig 6C, 6D and 6E. At the time point of 1200 s (outgrowing phase), the filopodium was found to start to protrude as its length was increased from 3.2 μm, the magnitude of elastic force ECM fibers, where filopodial tip are adhered, was found to be lowest peak as 72 pN (Fig 6A and 6B). It should be noted that simulated data were sampled at every ten seconds. In addition, lowest displacement fields of ECM fibers were observed at the local ECM area as the filopodia generate weak traction force through its outgrowing motion. Interestingly, directions of displacement fields are found to be away from the filopodia, which means ECM fiber network are relaxed.
At the time point of 1280 s (tugging phase), the filopodia was appeared to tug neighboring ECM fibers as its length was growing from 3.37 μm to 3.53 μm, the magnitude of elastic force ECM fibers, where filopodial tip are adhered, was found to be a high peak as 1.49 nN and the number of receptor-ligand bonds at filopodial tip was also found to be high as 248. In addition, highest displacement fields and substantial deformations of ECM fibers were observed at the local ECM area as the filopodia generate traction force through its tugging motion (Fig 6A and  6B and S12 Fig). Interestingly, directions of displacement fields are found to be towards the filopodia, which means ECM fiber network are tensioned.
At the time point of 1360 s (contractile phase), the filopodia was appeared to rapidly contract from 3.9 μm to 3.3 μm, the magnitude of elastic force ECM fibers, where filopodial tip are adhered, was found to be highest peak as 2.62 nN. In addition, highest displacement fields and substantial deformations of ECM fibers were observed at the local ECM area as the filopodia generate traction force through its contractile motion (Fig 6A and 6B). Interestingly, directions of displacement fields are also found to be towards the filopodia, which means ECM fiber network are tensioned.
At the time point of 1600 s (retractile phase), the filopodium was found to retract as its length was decreased to 3.2 μm, but the magnitude of elastic force ECM fibers was appeared to be a low value of 0.86 nN, and the number of receptor-ligand bonds was also moderately decreased to 215. Distinct behaviour of filopodial retractile motion is found to be the relaxation of tensioned ECM fiber network as the filopodia switches its motion from the contractile phase to retractile phase, which resulted in directions of displacement fields are observed to be away from the filopodium (Fig 6A and 6B and S12 Fig).

Discussion
It has been reported that distinct load-fail behaviour and frictional slippage behaviour of filopodia at both soft and hard gels can be explained by using the motor-clutch mechanism at the filopodial shaft [4]. Furthermore, the speed of the filopodia retraction is highly dependent on the stiffness of ECM surrounding the filopodia tip: as the stiffness of the ECM increases, the retrograded flow of actin at the filopodial shaft becomes faster, but the traction stress generated in the ECM becomes lower [4]. In the current stochastic model of filopodia penetration dynamics, these filopodia behaviours have been reproduced and their mechanisms have been elucidated as a dynamic 3D interplay between filopodia traction and ECM remodelling (S5-S7 Videos). Filopodia mechanically interact with ECM fibers to cause gel compaction and fiber remodelling [24] (S19 Video). In turn, the local stiffness of the ECM, which is modulated by the filopodia activities, influences the contractile/retractile behaviours of the filopodia via actin motor activities. Here, the contractile motion differs from the retractile motion in that the contractile motion leads to the generation of traction fields in the ECM surrounding the filopodia via the motor-clutch mechanism, while the retractile motion results in the relaxation of filopodia as well as the relaxation of ECM fibers. Furthermore, these distinct behaviours of filopodial motion have been confirmed using a 3D cell migration microfluidic assay in which 200 nmsized fluorescent beads are embedded into the collagen gel. As a result, beads move towards the extending filopodium tip during the tugging and contractile phases, and away from the filopodium during retractile motion (S13 and S14 Videos).
The ECM stiffness is known to increase as the concentration of crosslinking molecules increases [49]. On the other hand, an increase in the number of crosslinks leads to a decrease in the ECM pore size. In our current study, three kinds of ECM fiber network models were built with pore sizes of 0.5, 1.0 and 1.5 μm and single fiber diameters of 28, 34 and 41 nm, respectively, based on mechanical properties experiments [43]. The line stiffness of a single fiber ð¼ E e f A f =L e0 f Þ with a pore size for the three different ECM models were calculated to be 0.880, 0.908 and 1.2 pN/nm, respectively. Although these values are within the range of soft substrate whose stiffness is less than 1 pN/nm, apparent values of stiffness for the fiber network at the local area of 1 μm 2 are significantly increased to the range of harder substrate whose stiffness is more than 10 pN/nm as pore size becomes smaller and more crosslinks are added. To promote cell migration in a 3D ECM fiber network with the smallest pore size, degradation of the ECM fiber network is required [50]. For comparison we simulated filopodial penetration models with no degradation of the ECM fiber network, that is, MT1-MMP deficient cell model (S16 Video). This resulted in inhibition of the cell invasion into the ECM fiber network, compared to that observed for deep cell invasion model that incorporates filopodia penetration dynamics and ECM fiber network degradation (S7-S10 Figs and S2 Text and S17 Video). Thus, our simulated results reveal that the degradation of ECM fiber network plays an important role in filopodia penetration dynamics during both tugging and contractile phases. Without ECM degradation, filopodia can still grab ECM fibers; however, the lamella, where the root of the filopodial shaft is connected to the cellular membrane, can hardly penetrate through the ECM fiber network. With the addition of the component of degradation by MMP-2, our ECM fiber network model is able to simulate more complex situations; local fiber network, where MMP-2 is diffused from the cellular membrane, gradually becomes less stiffer since degraded fiber network is decomposed into multiple elements of single fibers whose stiffness are less than 1 pN/nm. Thereby, the lamella can easily penetrate through the degraded ECM fiber network.
Filopodia dynamics is very complex, and many other factors are likely to contribute to these dynamics. We have applied the Bell model to our filopodial focal complex model to rupture bonds between integrins and collagen molecules with force, and it is known as the 'slip bonds' as the ligand should slip out of the binding pocket more rapidly under higher tensile force [35]. As results, probability of rupturing an individual bond becomes higher and the lifetime of its bond becomes shorter under higher tensile force~200 pN. On the other hand, some cell-ECM interactions have been recently known to indicate 'catch bonds' behavior [51,52]; the lifetime of catch bonds show biphasic distributions under applied force, which implies that the lifetime of catch bond takes a maximum at a critical value of applied force. It is of interest how catch bonds likely to influence filopodia dynamics. There was an interesting work of stochastic models comparing behaviors between 'catch bond' and 'slip bond' in relation to actin retrograde flow and corresponding traction stresses [53]. In their work, in case of 'slip bond', bond fraction was rapidly decreased as actin retrograde speed was increased. However, in case of 'catch bond', bond fraction takes a maximum at actin retrograde speed of 5 nm/s. Overall, 'catch bond' case showed higher bond fraction than 'slip bond' case at the range of actin retrograde speed > 5 nm/s. Based on these results, filopodia dynamics with 'catch bonds' is likely influence longer lifetime focal complex formations than that with 'slip bonds'. Furthermore, filopodia dynamics with 'catch bonds' is likely influence stronger traction fields in ECM than that with 'slip bonds'.
The initial goal of current model was aimed at understanding how filopodia penetration dynamics plays an important role in 3D cell invasion into an ECM fiber network. The ultimate goal of the future cell model will be the development of the model like a realistic cell by capturing complicated morphology changes the cell. However, comparisons of cellular morphologies between the simulated cell models and experimentally measured HUVECs were limited in current cell model because morphologies of simulated cells were rounded or oval shapes, but those of experimentally measured cells were very elongated. To our knowledge, this discrepancy was resulted from two assumptions of the model; 1) the number of nodes in the cellular membrane was set to be 549, and 2) maximum number of clustered integrins per a node was set to be 100. Therefore, the direction of the future cell migration model will be to improve the morphology of cell model by increasing the number of nodes > 10,000 as the size of a mesh size is close to molecular size of~80 nm, and decreasing the number of clustered integrins. In addition, all diffusion coefficients and some secretion rates of biochemical concentrations in the current model were assumed to be identical for three ECM fiber models. However, in fact, heterogeneous diffusion coefficients and secretion rates of biochemical concentrations should be considered in the future cell migration model since recent experimental observations have indicated that diffusion coefficient of bovine serum albumin (BSA) is decreased as the ECM is stiffer and a pore sizes of ECM network is reduced more [22]. Furthermore, the future cell model can be extended to more complex models of cancer metastasis including a collective migration mediated by both cell-cell and cell-ECM adhesions, epithelial-mesenchymal transition (EMT)-mediated mesenchymal cell migration, amoeboid migration in ECM [54].

Intracellular mechanics
The structural mechanics and intracellular mechanics are other key mechanisms involved in cell invading into 3D gel. Here, we formulate them by extending the previous dynamic model for cell migration on curved surfaces [26], and cell migration and spreading on planar micropatterns [27]. The essential equations in the model include: 1) an equation for FA dynamics based on Monte-Carlo simulations of ligand-receptor bonds, 2) two equations for deformations of double elastic membranes: an outer cell membrane and an inner nuclear membrane, 3) an equation describing the contractile motion of actin stress fibers, which is extended from FAs on the cortical surface to the nuclear membrane, and 4) lamellipodium protrusion by actin polymerization [19] with a constant force of 300 pN. The detailed description of the equations in the model can be found in previous works [26,27]. Among them, the major extension in the mechanobiological dynamics model presented here is FAs dynamics in 3-D ECM fiber network model (Fig 1B). The FA force acts between the i-th integrin node on the cellular membrane and points of ECM fibers where the extension of the unit vector normal to the cellular membrane interacts with the nearest point of ECM fibers.

Construction of ECM fiber network model
To generate the computational model for ECM fiber network, a 3D cubic structure (50 × 50 × 50 μm) was initially drawn using AUTOCAD software. This structure was made in a stereolithography (STL) file format. Then, the surface geometry of the STL file was imported into a commercial CFD software package (STAR CCM++, CD-adapco) to build tetrahedral grids for the model. As a preprocess, further refinements to the STL file, including surface triangulation, tetrahedral volume meshing, and optimization for computational stability, were carried out to build three different ECM fiber network models with pore sizes of 0.5, 1, and 1.5 μm (S11A Fig). These preprocessed tetrahedral grids for three different ECM fiber network models with pore sizes of 0.5, 1, and 1.5 μm consisted of~3.0×10 6 , 3.8×10 5 , and 1.3×10 5 elements, respectively. Finally, these grids were modified to construct components of ECM fiber network, such as elastic fibers and crosslinks (S11B Fig). Here, red and yellow nodes in these grids represent crosslink and fiber nodes, respectively. To construct the fiber network, each line between two red nodes was divided by N div fiber segments and two crosslink segments that consist of N div +1 fiber nodes (yellow). Here, values of N div for ECM fiber network models with pore sizes of 0.5, 1, and 1.5 μm were respectively set to be 1, 2, and 3 with an assumption that the density of collagen molecules along all fibers in three ECM fiber network models was identical (300 collagen molecules per a fiber node), and the distance of a crosslink segment was set to be 30 nm. After segmentations of fibers, one set of fiber segments were randomly made to connect to coaxial neighbouring sets of fiber segments under a condition when the angle between two connected sets of fiber segments was above 60°(S11C Fig). In addition, the number of connected fibers at the i-th crosslink node was set to be P f , where P f is an initial ratio of forming fibers at the i-th crosslinks node (0.7), and N fs i is a total number of linked fiber segments at the i-th crosslinks node.

ECM fiber network dynamics
We assume the ECM fiber network to be composed of elastic ECM fibers and crosslinks, which make strong bonds between adjacent fibers [13]. The elastic energy stored in the ECM fiber network can be expressed in terms of the stretching and bending properties of the constituent fibers. The stretching modulus of a fiber is given by k e f ;s ð¼ E e f A f Þ, where E e f and A f ð¼ pr 2 f Þ are the Young's modulus (1 MPa) and the cross-sectional area of a single fiber, respectively. The bending modulus of a fiber is given by k e f ;b ð¼ E f I f Þ, where I f ð¼ pr 4 f =4Þ [55]. The stretching elastic energy of the j-th segment of the i-th fiber is given as a function of the difference between the stressed (L e ij ) and unstressed (L e0 ij ) lengths, and the bending elastic energy as the one of stressed (y e ij ) and unstressed (y e0 ij ) angles at the j-th node between two segments in the i-th fiber (Fig 7). The total elastic energy in the i-th ECM fiber in the network can be expressed as following: Here, it should be noted that the elastic energy at the j-th node in the i-th fiber is summed only for coaxial neighbouring nodes. Similarly, the elastic force at the j-th node in the i-th fiber,F e E;ij , can be derived by using the virtual work theory: where y e ik ¼ cos À1 ðt i;k Át i;kþ1 Þ;t i;k andt i;kþ1 are tangential unit vectors at the k and k+1-st segments in the i-th fiber, respectively, and To incorporate the degradability factor into the ECM fiber network and its nonlinear behaviour under mechanical responses, we consider that each crosslink node comprises crosslink molecules, such as amino acids, that can rupture. Fig 7B shows an example of connectivity between the i-th crosslink and the two neighbouring fibers. We model the degradability of ECM fiber network by considering detachment events among the i-th crosslink node and its neighbouring fibers, and the degradability of ECM fiber network depends on a local value of the ECM integrity (I ECM;i ¼ C ECM;i =C 0 ECM;i ; 0 I ECM;i 1) (see S18 Video). Here, C ECM,i and C 0 ECM;i are concentrations of the i-th ECM node at present and initial states (10μM), respectively. The number of uncrosslinked (or degraded) fibers at the i-th crosslink node,N uf i , is calculated as is an initial number of crosslinked ECM fibers at the i-th crosslinks node.

Numerical methods of "cell migration model in 3D ECM"
Cell migration simulations were carried out using a fourth order Rosenbrock method based on an adaptive time-stepping technique for integrating ordinary differential equations with the convergence criterion <10 −4 . The ordinary differential equations of cell and filopodia models were numerically coupled to solve for unknown variables associated with the mesh node position vectors for both cell membrane, nucleus membrane, transduce layer, and filopodial membrane (see Table 1). In particular, cell membrane and transduce layer were coupled with the viscoelastic actin cortex using kelvin-voigt model (Fig 1B). To solve two coupled ordinary differential equations in the CC and CT modules (Table 1) For cell migration simulation the Rosenbrock method outperforms the standard Runge-Kutta method which requires a relatively large number of iterations. Furthermore, the Rosenbrock method consumes less computing time by using adaptive time-step control that ranges from 10 −3 s to 10 −2 s in the present work. Thus, it is suitable for simulating transient cell migratory behaviours over 1 hour. For computations of ECM fiber networks, numbers of ECM fiber nodes were ranged from 30,000 to 234,000 depending on pore sizes of ECM model. As pore sizes of ECM fiber network models are smaller, computations of these models become more expensive. To solve ECM fiber network models effectively, computational domains of ECM fiber models were set within a radius of 15 μm at both the centre of cellular membrane and filopodial tip. As results, three ECM fiber network models (see Fig 4A, 4B and 4C) were visualized as half spheres. These computational domains were updated every 10 seconds of physical time as the cell interacts with ECM fibers dynamically.
Àk decay VEGF C VEGF MMP-2 Àk on TIMP2:MMP2 C TIMP2 C MMP2 þ k on Complex:MT1ÀMMP C Complex C MT1ÀMMP À k decay MMP2 C MMP TIMP-2 Àk on TIMP2:MMP2 C TIMP2 C MMP2 À k on TIMP2: Àk decay Ligand C Ligand þ k deg ECM C MMP2 C ECM X 2i R X2 i from the equation: One practical issue in computing finite mesh geometric models is to check geometrical compatibility. As the coordinates of cell membrane, filopodial membrane, and ECM fibril nodes are updated based on the equations of motion, geometrically incompatible situations occur occasionally in the configurations of the cell membrane mesh and that of the filopodia in relation to the ECM fibril surface. For example, some cell membrane nodes intersect with the fibril fiber, and the filopodia also intersects with the fibril fiber. These incompatible situations must be checked in every computational cycle, and necessary corrections such as contact forces (elastic repulsive) must be made [26,57].

Microfluidic experiments for filopodia dynamics
The methods for microfluidic sprouting experiments were previously described in detail [58]. Briefly, GFP-expressing HUVEC (Angio-Proteomie) were cultured in EGM2-MV growth medium (Lonza) and used in experiments at passage 6. PDMS microfluidic devices were bonded to glass coverslips and channels were coated with PDL (Sigma). The central gel region was filled with 2.5 mg/ml Collagen I solution before incubating for 30 min at 37°C to polymerize. The NaOH concentration of the Collagen I solution was varied to control the pH of the solution during polymerization. After polymerization of the gel, HUVECs were seeded in the medium channel by introducing 40 μl of cell suspension at 2×10 6 cells/ml. After 24 h, medium was replaced with EGM2-MV supplemented with 50 ng/ml VEGF (peprotech) and 250 nM S1P (Sigma). After 12 h, individual cells sprouting into the collagen gel were imaged at 60x magnification using a confocal laser-scanning microscope (FV-1000, Olympus). Collagen fibers were simultaneously visualized on the same instrument by collecting the reflected light (confocal reflectance microscopy). During imaging, cells were kept humidified at 37°C and 5% CO 2 . For imaging analysis, speeds of both filopodial tip and root was measured using a 3D image analysis tool (IMARIS, Bitplane). Here, L f is the length of filopodia, L max is the maximum length of outgrowing filopodia (4.5 μm), and L min is the minimum length of retractile filopodia (2 μm). T filo is the elapsed time of filopodia penetration dynamics. The rupture test indicates the procedure to determine the completed breakage of bonds at the filopodial tip by the Bell's equation. The polarity angle means the angle between direction of polarization for the cell and the unit vector normal to the local membrane where filopodial root is located. S2 Video. Filopodia tugging phase. Experimental observation shows that the binding sites of FCs advance along ECM fibers and that the bound ECM fibers are pushed or pulled by the filopodial tip. Green color indicates filopodial tip and shaft, and while lines indicate collagen fibers. This movie was processed using single sliced images. For full collapsed images of filopodia tugging phase, see S12 Video. (WMV) S3 Video. Simulation of filopodia tugging phase. Example of a simulated filopodium penetration into ECM network model with pore size of 1.5 μm over 970 seconds. Cell and filopodial membranes are visualized with green. Blue lines represent ECM fibers with single fiber's diameter of 41nm, and single fiber's modulus of 1 MPa. Scale bar is 500 nm. It shows that the filopodial tip crawls along ECM fibers, tugs neighbouring fibers, and contracts depending on the binding strength and stiffness and pore size of the ECM. S18 Video. Degradability of ECM fiber network model. An example of the degradability of ECM fiber network model due to the secretion of MMP-2 at the root of filopodium. Cell and filopodial membranes are visualized with green. Black lines represent ECM fibers with single fiber's diameter of 41nm, and single fiber's modulus of 1 MPa. Yellow dots represent crosslink nodes. Two blue arrows indicate clusters of crosslink nodes which are decomposed into multiple elements of single fibers. (WMV) S19 Video. Mechanical interactions between filopodia and ECM fibers. An example of gel compaction and fiber remodelling by filopodia penetration dynamics in ECM fiber network model with a pore size of 0.5 μm. Two yellow circles indicate significant gel compaction and fiber remodelling, and red arrows indicated directions of gel compaction toward to the cell. (WMV)