A computational modeling of invadopodia protrusion into an extracellular matrix fiber network

Invadopodia are dynamic actin-rich membrane protrusions that have been implicated in cancer cell invasion and metastasis. In addition, invasiveness of cancer cells is strongly correlated with invadopodia formation, which are observed during extravasation and colonization of metastatic cancer cells at secondary sites. However, quantitative understanding of the interaction of invadopodia with extracellular matrix (ECM) is lacking, and how invadopodia protrusion speed is associated with the frequency of protrusion-retraction cycles remains unknown. Here, we present a computational framework for the characterization of invadopodia protrusions which allows two way interactions between intracellular branched actin network and ECM fibers network. We have applied this approach to predicting the invasiveness of cancer cells by computationally knocking out actin-crosslinking molecules, such as α-actinin, filamin and fascin. The resulting simulations reveal distinct invadopodia dynamics with cycles of protrusion and retraction. Specifically, we found that (1) increasing accumulation of MT1-MMP at tips of invadopodia as the duration of protrusive phase is increased, and (2) the movement of nucleus toward the leading edge of the cell becomes unstable as duration of the retractile phase (or myosin turnover time) is longer than 1 min.

Cancer metastasis is a key step in cancer progression, and it is composed of multiple processes: (1) local infiltration of cancer cells into the adjacent tissue, (2) transendothelial migration of cancer cells into blood or lymphatic vessels (intravasation), (3) escaping from the immune system and survival in the circulatory system, (4) attachment of cancer cells to the vessel wall and migration out of the vessel (extravasation), and (5) subsequent proliferation in the distant organs leading to the establishment of colonies at the secondary sites 1,2 . Among these complicated processes, both intravasation and extravasation require invasive cancer cell migration into the basement membrane (BM) that is composed of highly specialized extracellular matrix (ECM). Cancer cell invasion into BM is facilitated by the dynamic movement of actin-rich cellular protrusion called invadopodia 3 . Hence, examining how cancer cells utilize actin-myosin machinery during dynamic movement of invadopodia protrusion (IP), as well as how they preferentially change their morphologies during the extravasation, will provide valuable insights into the intricacies of cancer metastasis.
The formation of invadopodia is a characteristic of highly invasive cancer cells 4 . Invadopodia are elongated ventral membrane protrusions that are composed of a variety of proteins, such as actin filaments, actin-related protein-2/3 (Arp2/3) complex, neuronal Wiskott-Aldrich syndrome protein (N-WASP), cortactin, fascin, and matrix degradation enzymes. Among them, N-WASP and cortactin are essential components that can synergistically activate Arp2/3 complex, and these proteins are upregulated in malignant cancer cells 5 . Invadopodia play a critical role in the process of extravasation by inducing basement membrane (BM) disruption through local ECM degradation. In addition, it has been reported that invadopodia is potently induced by high-density fibrillary collagen (HDFC) matrix in both cancer cell lines and primary human fibroblasts 6 . To achieve this, cancer cells need to secrete matrix metalloproteinase (MMP) family to degrade local ECM with which cellular membrane interact. Specifically, matrix type 1-metalloproteinase (MT1-MMP) accumulates at tips of invadopodia, which facilitates membrane protrusion. On the other hand, an interesting study has recently revealed a new mode of invadopodia-assisted migration that is ECM plasticity-mediated and protease-independent. This

Results
Cancer cell invasion model with invadopodia. The activity of invadopodia, as well as its interaction with 3D ECM 13 , is a complicated, multifaceted process where at least three individual dynamics are coupled: (1) cell mechanics, including the protrusion and retraction of invadopodia membrane and the growth of discrete intracellular branched actin network; (2) mechanics of the discrete ECM fiber network (S1 Text); and (3) reaction-diffusion mass transfer over ECM (Supplementary Fig. 1 and Supplementary Information). In addition, to incorporate viscoelastic behaviors in (1) triple layers of cellular membrane (invadopodia membrane, force transduction layer (FTL), and actin cortex layer (ACL)), (2) intracellular branched actin network, and (3) double layers of nuclear membrane (perinuclear actin layer (PAL), and nuclear membrane surface (NMS)) ( Fig. 1a), all the line elements in their structures were simulated with Kelvin-Voight model (a spring and a dashpot in parallel). Figure 1b shows the free body diagram of the i-th node on the invadopodia membrane, called the i-th integrin node, where a bundle of integrins is formed. The detailed equations that govern each of these dynamical processes and a list of simulation parameters are summarized in Table 1 and Supplementary Table 1, respectively. Invadopodia dynamics. Dynamic behaviors of invadopodia and their interactions with the network of intracellular F-actins are complex, but can be predicted with a computational model based on a few assumptions and observations. First, the discrete branched actin network is composed of viscoelastic actin filaments (diameter of 7 nm and modulus of 1.8 GPa) 14 , Arp2/3 complexes, capping proteins (CP), actin crosslinking molecules (α-actinin, filamin, and fascin), and contractile bipolar myosin filaments (Fig. 2a-c). Second, the invadopodia protrusion dynamics into 3D ECM consist of three different motion phases in a single cycle (Fig. 2d) 12 : 1) a protrusive phase driven by F-actin polymerization, 2) a retractile phase due to contractile motions of perinuclear bipolar myosin filaments, and (3) a severing phase 15 for the rapid depolymerization of buckled actin filaments 16 by actin severing proteins (Cofilin and Gelsolin). Particularly, the duration of protrusive phase has a key role in facilitating invadopodia protrusion into ECM since MT1-MMP, which is accumulated at tips of invadopodia, can degrade surrounding ECM only during this phase. Accordingly, MT1-MMP signal pathway is assumed to be switched on only during this phase.
Protrusive phase in discrete F-actin mechanics. Here www.nature.com/scientificreports/ growing rate and v s is shrinkage rate of 10 nm/s 17 . v g is constant or variable depending on the adhesion of barbedend for actin filament to the actin cortex layer (ACL). v g is constant (20 nm/s) when a daughter actin filament is nucleated from a mother filament. Then, as growing actin filament starts interacting with the surface of ACL, v g becomes variable. To model the mechanics of growing actin filament against the surface of ACL, concave-up velocity-force relation (exponential decrease of the velocity with the growing load) is used 18 . Here where v 0 is 50 nm/s, f cortex is the tension at a node on the ACL, n F is the number of actin filaments interacting with the node on the ACL, and f 0 is 20 pN. Thereby, the effective stiffness of the j-th line elements of the i-th actin filament is variable, , where E A is the Young's modulus of actin filament (1.8 GPa) 14 and A A is the average cross-sectional area of actin filament (38.48 nm 2 ). It should be noted that an initial value of L A0 i,j is 150 nm, and thus an initial value of κ A L,i,j is 0.46 N/m. In our previous publication, we applied the virtual energy method in structure mechanics to derive elastic forces at the cellular and nuclear membranes 19,20 , and ECM fibers 21,22 . Here, in a similar manner, we apply the virtual energy method to derive elastic forces at actin network include F-actin, crosslinking molecules, and myosin minifilaments. In addition, a previous study used MEDYAN model and energy potential method to calculate elastic forces due to stretching and bending polymer as well as branching and dihedral angles 23 . It should be noted that our method is different from theirs in that we analytically derive elastic force by differentiating total elastic energy at the specified node. Here, the total elastic energy stored in the i-th actin filament is given by where L A i,j is the length of the j-th line of the i-th actin filament and is updated at every time-step are location vectors of the j-th and j + 1-th nodes of the i-th actin filament) (Fig. 3a). L A0 i,j is its relaxed length (150 nm). The other is the total elastic energy associated with the change of angle between L A i,j and L A i,j+1 on the i-th actin filament, and it is expressed as where κ A b is the bending modulus of actin filament is expressed as i ,j are the stressed and unstressed angles (zero) between L A i,j and L A i,j+1 on the i-th actin filament, respectively  www.nature.com/scientificreports/ ( Fig. 3a). Similarly, the elastic force at the j-th node on the i-th actin filament, F A E,i,j , can be derived by using the virtual work theory: , t i,k and t i,k+1 are tangential unit vectors at the k and k + 1-th lines in the i-th actin filament, respectively, and In case that a new k-th actin filament (daughter) is nucleated from the Arp2/3 complex which is on the j-th node on the i-th actin filament (mother), an elastic Arp2/3 force,F A arp,i,j , and elastic branch force,F A br,i,j , are nonzero (Fig. 3b). Similarly,F A arp,i,j is expressed as where κ arp is the effective stiffness constants of Arp2/3 complex (0.01 N/m), L A arp,i,j is the distance between the j-th node on the i-th actin filament and the first node on the k-th actin filament, and arp is unstressed length of   Table 1 can be found in Supplementary Information.

Module Equation Couplings
CI n(i) www.nature.com/scientificreports/  www.nature.com/scientificreports/ Arp2/3 complex (30 nm). F A br,i,j has two kinds of branch forces: one is branch angle force, and the other is branch dihedral force 23 . Total branch elastic energy at the j-th node on the i-th actin filament is expressed as where κ A br,ang is the angular bending modulus 23 , φ A i, j is the branched angle between L A i,j and L A k,1 . That is, the angle between two unit vectors of l i,j = , and φ A0 i ,j is an equilibrium branched angle ~ 70˚ between the mother and daughter filaments. κ A br,dihed is the dihedral bending modulus 23 , ψ A i, j is the dihedral angle between two planes, formed by the points , and ψ A0 i ,j is assumed to b e z e ro. Two ve c tors n or m a l to t h e pl an e s are re sp e c t ive ly e x pre ss e d a s ×l k,1 . Thereby, F A br,i,j can be derived by using the virtual work theory: where φ A i, j = cos −1 l i,j ·l k,1 , and ψ A i, j = cos −1 m i,j ·m k,1 . Similarly, Furthermore, a polymerization force due to Brownian ratchet motion 18,24 , which is induced by thermal fluctuations of membrane, is also incorporated in the model, and it is expressed as where F p is assumed to be 20 pN because the force exerted by each actin filament is the order of 10 ~ 20 pN 25 , v 0 is zero load polymerization speed of actin filament (50 nm/s), n I R,i is a unit vector normal to the surface of i-th www.nature.com/scientificreports/ invadopodial membrane, v b is a velocity vector at the barbed-end for actin filament (CA module), and v a is a velocity vector at the node of actin cortex layer mesh (CC module) where barbed-end for actin filament interacted with. It should be noted that this polymerization process due to Brownian ratchet motion is switched on during the protrusive phase, but it is switched off during both the retractile and severing phases since binding of capping proteins (CPs) at barbed-ends for actin filaments block the addition or loss of G-actins.

Retractile phase in discrete F-actin mechanics.
To incorporate the retractile phase in the invadopodia dynamics, mechanical interactions of 3D branched actin filaments by actin crosslinking molecules (fascin, filamin and α-actinin) and contractile bipolar myosin filaments are modelled. We assume that fascins connect two parallel adjacent actin filaments which are bundled in the core of invadopodia, and filamins connect two parallel or orthogonal actin filaments which are located in the leading edge of the cell. It should be noted that both fascins and filamins were assumed to form in the leading edge of the cell, where their directional angles between a unit vector of cellular polarity (migration direction) and their location vectors from the center of cell were less than 60 degrees. In case of α-actinin and bipolar myosin filaments, they are assumed to be activated during the retractile phase and connect two parallel actin filaments which surround the nucleus 26 23 . Here α ∈ [0, 1] and β ∈ [0, 1] represent stochastic fractional ratios. Thereby, linking elastic forces at the j-th node on the i-th actin filament,F A L,i,j , can be derived by using virtual energy method, and it is expressed as Similarly, the linking elastic force at the m-th node on the k-th actin filament, F A L,k,m , can be expressed as It has been known that a bipolar myosin filament move towards barbed-ends for two adjacent actin filaments. The elastic energy at the l-th bipolar myosin filament, which connect the j-th line element on the i-th actin filament and the m-th line element on the k-th actin filament, is expressed as where κ BMF C is the effective stiffness constant of the bipolar myosin filament (0.01 N/m), and L BMF and L BMF0 l are the stressed and unstressed lengths of the l-th bipolar myosin filament (300 nm) 27 , respectively. Similarly, x A i ,l and x A k ,l indicate binding positions on the j-th line element on the i-th actin filament and the m-th line element on the k-th actin filament, respectively. They are expressed as . Similarly, contractile force at the j-th node on the i-th actin filament by a bipolar myosin filament, F A C,i,j , can be expressed as Contractile force at the m-th node on the k-th actin filament by a bipolar myosin filament, F A C,k,m , can be expressed as It should be noted that α and β are updated at every time-step because of sliding motions along to two actin filaments. Let x A * i ,l be the predicted binding position at the next time, x A i ,l represents the biding position at the current time, and it can be expressed as where v sliding is the sliding rate of bipolar myosin filament, and △ t is the time-step. Let α * be the predicted α at the next time-step, and it is expressed as www.nature.com/scientificreports/ The above equation can be further expanded and simplified as following: In addition, to incorporate mechanical interplay between the non-muscle myosin II and actin filaments, we adopt the force-velocity relation of muscle myosin II, first proposed by A.V. Hill 28 , and v sliding in Eq. (12) can be expressed as where v sliding0 is the sliding rate of myosin in the absence of load (160 nm/s), F stall is the stall force (1.24 Nn), c m is a dimensionless myosin parameter (4.761) 32 , and F ext is the magnitude of sensed contractile force at the previous time-step. The retractile phase is related with the disassembly of non-muscle myosin II with alphaactinins. At the end of the retractile phase, α-actinins were assumed to turnover to terminate the contractile activity of non-muscle myosin II. Thereby, the time duration in the retractile phase means turnover time of α-actinins. In addition, we also modelled the turnover of fascin and filamin to incorporate dynamic behavior due to their rupturing and rebinding to two adjacent actin filaments. After parameter study on turnover time of α-actinin (see Fig. 4c), we assume to set turnover times of crosslinking proteins as 60 s, which is reasonable since measured lifetimes of the α-actinin/actin interactions range from 15 to 155 s. For example, by molecular rupture measurement under constant load, lifetime of the α-actinin/actin interaction was found to be ~ 20 s (corresponding to k actinin off ≈ 0.05s −1 ) 29 . Interestingly, intrinsic dissociation rates for α-actinin and filamin, which were estimated with theoretical model using overall rupture-force probability distributions, are 0.066 ± 0.028 s −1 and 0.087 ± 0.073 s −1 , respectively 30 . Furthermore, lifetime of the disease-causing mutant K255E of human α-actinin-4 was measured as (155.1 s) 31 .
Severing phase in discrete F-actin mechanics. Severing phase has a key role in the rapid depolymerization of buckled actin filaments 16 resulting from contractile motions of perinuclear bipolar myosin filaments. As actin severing proteins (Cofilin and Gelsolin) can enhance disassembly of single actin filament by creating multiple pointed-ends of segmented actin filaments, a parallel mode of depolymerization process occurs at these pointed-ends of multiple actin filaments simultaneously. Thereby, severing phase with the parallel mode of depolymerization is faster than a serial mode of depolymerization. Here, we computationally implement severing process by segmenting single actin filament to multiple actin filaments with a length of 300 nm. Then, the lengths of segmented actin filaments are computed with the equation, dL A0 i,1 dt = −v s at the pointed-end of actin filament where v s is shrinkage rate of 10 nm/s. We assume that severing duration time is set to be 20 s.
Dynamic state change of F-actin. During the dynamic process, the state of single actin filament changes depending on the adhesion of barbed-end for actin filament to the actin cortex layer (ACL). This can be modelled as a discrete state transition network detailed in Supplementary Fig. 2. Briefly, a daughter actin filament is nucleated at Arp2/3 complexes on a mother actin filament with a nucleated rate (k n ) of 50 s −1 (nucleated state, S = 1), then the daughter actin filament grows with a directional angle of 70 degrees from the mother actin filament (growth state, S = 2) until the barbed-end for actin filament contact to the surface of the actin cortex layer. In addition, the formation of adhesions of actin filaments to the surface of ACL is assumed to occur when h p (gap between the barbed-end for actin filament and the surface of ACL) is less than 100 nm (adhesion state, S = 3). This adhesion is also assumed to rupture when tensile force is higher than 20 pN (detach state, S = 4) due to Brownian ratchet motion at the cellular membrane, and then immediately actin polymerization occurs (ratchet state, S = 5). In particular, it is known that increasing the capping protein (CP) concentration decrease the density of actin network, while increasing the Arp2/3 complex concentration increases the network density 33 . Accordingly, during invadopodia retraction and severing phases, all the barbed-ends for actin filaments at the tip of invadopodia are forced to bind with CPs to block the actin polymerization (capping state, S = 6). For the rest of leading edge of the cell, CPs are assumed to bind to barbed-ends for actin filaments (growth or ratchet state) with a binding rate (k cap ) of 0.3 s −1 . Lastly, CPs are uncapped to resume actin polymerization when invadopodia protrusion phase is recycled (uncapping state, S = 7).
Characterization of invadopodia protrusion. First, three series of simulations were performed to investigate the effects of protrusive phase on the protrusion of invadopodia into ECM by changing time durations of protrusive phase (60, 120, and 240 s) (Fig. 4a,b). Here, time durations in retractile and severing phases for each case were kept constant at 60 and 20 s, respectively. The time duration in the retractile phase represents the turnover time of non-muscle myosin II. Computational model of branched actin network was constructed with two patterns of actin filaments. Initially, the first pattern of 500 actin filaments as a base case was randomly aligned so that its barbed-end and pointed-end for actin filaments are close to the inner cellular membrane and the nuclear membrane, respectively. Each actin filament can rapidly polymerize at the bared-end with a rate of (13a) www.nature.com/scientificreports/ 20 nm/s, depolymerize at the pointed-end with a rate of (10 nm/s) 23 , and Arp2/3 complex binds the mother filaments and initiates the polymerization of daughter filament at a branch angle of 70 degree from the mother filament. Then, the second pattern of 500 actin filaments was also randomly distributed to encircle the nuclear membrane, and a pair of these filaments were cross-linked by a bipolar myosin filament with an unstressed length of 300 nm and α-actinin with an unstressed length of 50 nm. Computational model of ECM fiber network (45 × 30 × 20 µm) is constructed with a pore size of 3.0 µm and single collagen type 1 fiber diameter of 41 nm. The presented computational method of ECM and its simulated results are consistent with recent cell migration experimental observations. Computational model successfully reproduced these observations that speeds of both tip and root of filopodia increase as the pore size the collagen gel is increased in experiments 21 . Previously, to characterize bulk modulus of ECM fiber network with a pore size of 3.0 µm, a simulation of mechanical stretching test for corresponding ECM fiber segment was performed, and bulk modulus of ECM model was found to be (2,558 Pa) 22 . Then, spherical cell model was initially placed in the center of ECM domain, and two kinds of volume exclusion effects of the cell were considered (i) to prevent ECM fibers from penetrating though external cellular membrane, and (ii) to keep F-actin in the intracellular domain with two surface boundaries of internal cellular membrane and nuclear membranes ( Supplementary Fig. 3). As a result, bipolar myosin filaments are densified and aligned along the surface of the nuclear membrane. Thereby, the nucleus was pulled toward to the branched actin filament network because strong contractile force www.nature.com/scientificreports/ was generated by densified bipolar myosin filaments ( Fig. 4a and Supplementary Movie 1). In addition, we found that as the duration time in the protrusive phase increases, the formation of invadopodia protrusion was enhanced since MT1-MMP can be secreted to degrade ECM surrounding the cells (Fig. 4b and Supplementary Movie 2). We also found that the speed of nucleus was faster when the time duration in the retractile phase (or turnover time of non-muscle myosin II) is longer (Fig. 4c). Interestingly, multiple cycles of invadopodia elongation and retraction is captured for all three cases (Fig. 4d), reproducing the experimental observation 12 . The simulated migrated distance of invadopodial tip showed an excellent linear correlation to the duration time of protrusive phase (Fig. 4e). Second, three series of simulations were performed and compared to investigate the effect of initial number of F-actin (1000, 2000, and 3000) on the invadopodia protrusion ( Supplementary Fig. 4). In these simulations, time durations in the protrusive, retractile, and severing phases were kept constant at 240, 60, and 20 s. At these parameters, the speed of invadopodia was highest according to the first set of simulations. At time-point of 900 s, simulated plots of all three cases shows sharp invadopodia protrusions with densified branched actin network into ECM fiber network (Supplementary Fig. 4a-c). Hence, both the number of adhered F-actin to the cellular membrane and the speed of invadopodium increase as the number of initial F-actin increases ( Supplementary  Fig. 4d,e). Statistical analysis of linear regression was performed by comparing the number of adhered F-actin and the speed of invadopodium in terms of the number of initial F-actin. Excellent correlations were found between the two with r 2 = 0.998 (Supplementary Fig. 4f.).

Experimental observation of invadopodia protrusion and MT1-MMP inhibition.. We elect to
test the validity of our computational model by observing the dynamics of invadopodia formation of cancer cells cultured within a 3D collagen I extracellular matrix (ECM). We chose MDA-MB-231 breast carcinoma cells in our study as they were widely reported to produce invadopodia when cultured in ECM [34][35][36] . To observe the interaction of the invadopodia with the collagen type I ECM, we perform live-cell time-lapse confocal microscopy imaging of cancer cells cultured in the microfluidic device. MDA-MB-231 cancer cells cultured in collagen type I ECM are elongated and produce invadopodia to probe the surrounding matrix. The time-lapse images of the invadopodia reveal that they are incredibly active, and cancer cells use invadopodia to probe the surrounding matrix (Supplementary Movie 3). To quantify the invadopodia protrusion dynamics, the movement of invadopodial tip from each cell was tracked, and the average speed of this movement was calculated to be 22.5 nm/s.
To evaluate our in-silico model, we incubated MDA-MB-231 cancer cells with blocking antibody against MT1-MMP, and quantified the invadopodia movement speed as described before. Inhibition of MT1-MMP activities by blocking antibody significantly reduced the activity of the invadopodia (red arrow in Fig. 5a, and Supplementary Movie 4) and the average invadopodia speed (from 22.5 nm/s to 10.5 nm/s, Fig. 5b). Moreover, the distribution of the invadopodia speeds shifts toward lower values as the result of MT1-MMP inhibition (Fig. 5c). These results match well with the outcome of our computational model with MT1-MMP knockout and the inhibition of Arp2/3 complexes (Fig. 5d and Supplementary Movie 5), indicating more significant reduction of the average invadopodia speed at this model (13.4 nm/s to 6.8 nm/s, Fig. 5e) than that of a model with MT1-MMP knockout only.

Cancer cell invasive models with knockout of α-actinin, and filamin and fascin.. It has been
known that overexpression of filamin A has a tumor-promoting effect when it is localized to the cytoplasm 37 . In addition, fascin up-regulation in the invadopodia is highly related to increased cancer cell motility 38 , and increased α-actinin-1 level promotes cell migration, especially in human breast cancers 39 . On the other hand, fascin elimination decreases migration in glioma cells 40 and downregulations of both α-actinin-1 and α-actinin-4 make critical and distinct contributions to cytoskeletal organization, rigidity-sensing, and motility of glioma cells 41 . Therefore, cancer cells with impaired α-actinin, filamin, and fascin may have a reduction in motility. To investigate how α-actinin, filamin, and fascin affect cancer cell invasion (Fig. 6a,b), we computationally knocked out these proteins in our simulation and observed the resulting effects on invadopodia protrusion into ECM. Myosin filaments were observed to scatter out of the nucleus (Fig. 6a), which resulted in weak contractile force transmitted to both invadopodial and nuclear membranes (Supplementary Movie 6). Interestingly, invadopodia formation was significantly inhibited (Fig. 6b and Supplementary Movie 7), which resulted from the unbundling of actin filaments. Hence, knockouts of filamin and fascin significantly reduced the activity of the invadopodia (Fig. 6c,d) and the averaged invadopodia speed.

Time evolutions of extracellular and intracellular forces during the directed cell migration..
Comparing the series of simulated data of actin filaments and invadopodia provides insights into the time evolution of traction and intracellular forces during the directed cell migration toward stiffer ECM. To analyze dynamic behavior of invadopodium during stiffness-directed cell migration, the simulation was set up with cells located 10 µm from the sharp interface between soft and stiff ECMs. In this case, time durations of protrusive, retractile, and severing phases were set to be 300, 60, and 20 s, respectively. It has recently been reported that the secretion rate of MT1-MMP at the tip of invadopodium is enhanced in a stiffer ECM 13,42 . Accordingly, in this simulation, the secretion rate of MT1-MMP at the tip of invadopodium is assumed to be positively correlated with the density of surrounding ligand. As shown in Fig. 7a,c, three selected plots of branched actin network and ECM fiber network at time points of 300 s, 355 s, and 560 s show how branched actin filaments are gradually densified as the tip of invadopodium approaches stiffer ECM (Supplementary Movie 8). Temporal variations of traction force (Fig. 7d), intracellular force (Fig. 7e), and trajectory (Fig. 7f) at the tip of invadopodium were evaluated in the time range of 100-500 s. Time-averaged traction force and intracellular force were calculated as 0 .10 nN and 1.16 nN, respectively. In Fig. 7f,  www.nature.com/scientificreports/ podia elongation and retraction, with maximum at 260 s and minimum at 380 s. Unlike the first cycle, protrusive distance of invadopodial tip at the second cycle was reduced because of sharp interface between soft and stiff ECMs, which indicates significant increase in traction force at the tip of invadopodium. In this report, our simulated data reveal that contractile forces, adjusted by the turnover time of perinuclear myosin filaments, can induce the movement of nucleus, which can facilitate the cancer cell invasion and migration into extracellular matrix confinement. Interestingly, a recent study revealed that non-muscle myosin II filaments have fast turnover with a half time of 1 min, while the half-life in skeletal muscle cells is more than an hour 43 . Accordingly, we reflect this fast turnover time into our model, and found that the fast turnover induces nuclear movement towards the branched actin network. As for the nuclear movement, it has been reported that the formation of perinuclear actin fibers induces temporal rotational movement of the nucleus, which results in nuclear reorientation to the direction of migration 44 . Further, our computational model was able to predict the role of actin crosslinkers in invadopodia formation by computationally depleting α-actinin, filamin, and fascin, from the branched actin network. For example, through a simulation of α-actinin depletion, we reveal that bipolar myosin filaments are impossible to generate contractile forces without α-actinin, resulting in impaired nuclear movement. In addition, in case of a simulation for the knockout of filamin and fascin, we also found that invadopodial membrane cannot protrude without bundling of actin filaments by filamin and fascin.
We introduce a computational framework for simulating the invasion of metastatic cancer cell into ECM by integrating mechanics of intracellular branched actin filaments and discrete ECM fibers in the invadopodia protrusion dynamics. This computational framework includes two-way interactions between intracellular and extracellular domains. The effects of intracellular domain on extracellular one are represented by the remodeling and densification of ECM via mechanical force generated by invadopodia and chemical degradation of (where p is perimeter of single cell, A is the area of the cell), is larger than 3.81, which represents fluid-like tissue or unjammed cells of disordered metastable tissue configurations 45 . For example, the cell with many filopodia behaves like fluid since its perimeter is larger. However, in our case, we assume that invadopodia protrusion has increased actin branching density than that of filopodia protrusion 46 . MT1-MMP accumulates at tips of invadopodia and plays a crucial role in focal pericellular degradation of the BM 47 . Furthermore, it has been proposed that high invasiveness of cancer cell migration occurs at the intermediated levels of ECM crosslinking (0.39) which strictly depends on proteolysis of the matrix by MT1-MMP 48 . It is known that the activation of integrins leads to the assembly of F-actin, Arp2/3 complex, formins, and cortactin in the invadopodia 49 . Aggregation of F-actin and cortactin initiates the accumulation of MT1-MMP at the invadopodia. F-actin density and MT1-MMP are related because experimental findings show that depletion of the protease impairs F-actin and cortactin accumulation at the invadopodial tip 50 . In particular, cortactin is one of key components involved in the delivery of MT1-MMP to invadopodia, and cortactin has an important role as a regulator of arp2/3-mediated actin branching 51 . Furthermore, Arp2/3 overexpression is tightly associated with cancer progression and tumor cell invasion. Arp2/3 overexpression was also associated with poor patient survival in lung 52 and breast cancers 53 . Arp2/3 overexpression in breast cancer is associated with HER2 overexpression 54 . Based on above experimental findings, our simulations suggest that Arp2/3 is involved with the secretion of MT1-MMP, which is not necessarily verified. Accordingly, we reflect this correlation between Arp2/3 mediated actin branching and the secretion of MT1-MMP in our model to match our computational www.nature.com/scientificreports/ simulation with experimental results of MT1-MMP knockout. Thereby, we found more significant reduction in average invadopodia speed using our current model, which considers both MT1-MMP chemical signaling and Arp2/3-mediated mechanical branching of actin (Fig. 5e), than a model that only considers MT1-MMP signaling. Current model for the invadopodia protrusion is inspired by above-mentioned experimental findings. We assume that single F-actin can induce Arp2/3-mediated actin branching when this F-actin is connected to the activated integrins through the force transduction layer. Subsequently, actin branching is further densified at the proximal invadopodium, which results in the recruitment of more activated integrin clustering and strong traction force and stress over the surrounding ECM. Recently, an interesting finding is the adaptation of branched actin network, that is, actin branching become denser when higher external force is exerted on actin network 55 . Thereby, as the surrounding ECM is stiffer and denser, the size of focal adhesions is larger, which will result in denser branched actin network. Eventually, cell needs to modulate to high secretion rate of MT1-MMP. In this aspect, the two-way interactions between cell (actin) and ECM can be coupled and feed-backed each other to promote invadopodia activity.
We note that our model of the invadopodia protrusion into ECM fiber network is complex and multifaceted. This model takes into account of impacts ECM stiffness and density on branched actin network, initial concentrations of actin and crosslinking molecules, rates of actin polymerization and depolymerization at the barbed-ends www.nature.com/scientificreports/ and pointed-ends for actin filaments, turnover time of mechanosensitive contractile bipolar myosin filament, and binding rate of CP at the barbed-end for actin filaments. Among them, the most sensitive parameter is the turnover time of myosin since turnover time > 90 s results in strong contractile force that pulls the nucleus too close to the invadopodial membrane. A recent study found that a rigid ECM promotes MT1-MMP secretion at the invadopodia 13 . Importantly, both the adaptation of branched actin network and MT1-MMP secretion rate are controlled by the surrounding ECM stiffness, which represents invadopodia mechanosensing. Current computational model is limited in feedback control of MT1-MMP secretion by cell-ECM interactions. Therefore, future invadopodia model will need to incorporate this invadopodia mechanosensing. In addition, current model is also limited in the characterization and validation of single actin network model only, for example, the dependence of viscoelastic properties of actin network on crosslinker concentration and myosin activity. Instead, we looked into whether dynamic behaviors of invadopodia depend on myosin activity or actin crosslinking molecules by depleting these molecules in our simulations. To our knowledge, this is the first 3D model that predicts the experimentally observed invadopodia behaviors in respond to the inhibition of actin-crosslinking molecules.
In addition, we focused on the validation of whole comprehensive model of cell-actin-ECM because it would be better to provide an insight of the model. Therefore, the extension of our current model will be to simulate how actin network will respond to indentation by deflectable cantilever and to predict the experimental results obtained via atomic force microscopy measurements 56,57 . Our current membrane model has fixed topologies and has a limitation in forming very narrow protrusions like filopodia. To overcome this limitation, we previously constructed filopodia geometric model by allowing the growth of filopodium patches with six triangular elements (or a hexagonal patch) on the cellular membrane 21,22 . These triangular elements were remeshed and new elements were allowed to generate along to the filopodium shaft. Interestingly, there has been a novel grand canonical membrane simulation model that describes stochastic polymerization of filaments against a fluctuating fluid membrane 58 . Their results suggest a membrane mediated dynamical transition from single filaments to cooperatively growing bundles as an important dynamical bottleneck to tubular protrusion. Finally, our ultimate goal is to use the computational model to understand the role that invadopodia protrusions play in cancer cell invasion and metastasis. To achieve this, the future direction of this model will need to simulate cancer cell extravasation and subsequent invasion into thin but dense BM along the blood vessels.
In summary, there have been few comprehensive models of invadopodia protrusion and invadopodia-ECM interactions due to the complexity of both the branched actin network and ECM fiber network. Moreover, how invadopodia speeds are related to the frequency of invadopodia protrusion is also poorly understood. Using the computational model we proposed in this study, we found that invadopodia speeds become faster as the protrusive phase becomes longer. In addition, using this computational model to simulate the knockout of actin-crosslinking molecules, we evaluated the effects of α-actinin on perinuclear bipolar myosin filaments, and filamin and fascin on branched actin network. Specifically, we showed that 1) it is impossible to generate contractile force without α-actinin, resulting in impaired nuclear movement; and 2) bundling of actin filaments by filamin and fascin is required for the invadopodia protrusion. The computational model was then applied to test why non-muscle myosin has a fast turnover time. We found that myosin turnover time of 1 min is enough to produce stable movement of the nucleus towards the leading edge of the cell. Altogether, we believe that our computational model has the potentials to be a critical tool in understanding cancer invasion and metastasis in various ECM microenvironments.

Methods
Computational model of invadopodia protrusion dynamics. The elastic force at the i-th node of invadopodia membrane (module CI in Table 1), F I E,i , is derived by using the virtual work theory in structural mechanics. To this end, the total elastic energy stored in the invadopodial membrane is obtained. Two types of total elastic energies are considered. One is the total elastic energy associated with distance changes between the nodes 59,60 : where L I i is the length of the i-th line of the invadopodial membrane mesh and is updated at every time-step. L I0 i is its relaxed length. κ I L is effective stiffness constants of the line elements of the invadopodial membrane (5.0 × 10 -5 N/m) 61 . The other is the total elastic energy associated with area changes: where A I i is the i-th mesh area of the invadopodial membrane and A I0 i is its relaxed values. κ I A is an effective stiffness constant of area elements of the invadopodial membrane (1.0 × 10 -4 N/m 2 ) 61 . Then, F I E,i can be obtained by differentiating these two kinds of total elastic energy, The focal complex force at the i-th node of invadopodia membrane,F I FC,i ,is expressed as www.nature.com/scientificreports/ where n I b,i is the number of integrin-collagen bonds at the i-th node of invadopodial membrane,κ LR is the spring constant of a single integrin-collagen bond (~ 1 pN/nm) 62 , L b is the average stretched length of the integrincollagen bonds, is an unstressed length of bonds (~ 30 nm) 63 and n I R,i is a unit vector at the local surface of the i-th node of invadopodial membrane toward the bonding site at the ECM fiber. Here (L b − ) represents the stretched distance from the equilibrium. We utilize Bell's model 64 to incorporate force-dependent reaction rates of n I b,i is expressed in following ordinary differential equation: where n I tot is total available number of integrin molecules at the i-th node of invadopodial membrane, k on is the kinetic associate rate for binding integrin molecules and ECM fiber, and it is expressed as 65,66 where k 0 on is the zero forward reaction rate (1 molecule −1 s −1 ), l bind is a binding radius (30 nm) to check whether the i-th node of invadopodial membrane and the node on the fiber are sufficiently close, and k b T is the unit of thermal energy.Z 0 is the partition function for a integrin molecule confined in a harmonic potential between − and L b − , and it is expressed as k off is the kinetic dissociation rate, and it is known as Bell's equation for the slip bond, which is defined by 64 where k 0 off is the zero kinetic dissociation rate in the absent of the force, x b is the distance between the minimum binding potential and the transition state barrier, and k b T x b represents an intrinsic force ~ 200 pN. The transduce force at the i-th node of invadopodia membrane,F I T,i , is expressed as where κ T is an effective spring constant of line element of the FTL (8.0 × 10 -3 N/m), L T,i is the tensioned length of the i-th line in the FTL, and it is updated at every time-step. L 0 T,i is its relaxed (zero force) length (50 nm). The elastic force at the i-th node of force transduce layer (FTL) (module CT in Table 1) The elastic force at the i-th node of ACL (module CC in Table 1),F ACL E,i , is expressed as where A A is the sectional area of the actin filament ( πr 2 A ), where r A is the radius of actin filament (3.5 nm), and E A is Young's modulus of the actin filament, L C i is the length of the i-th line on the ACL mesh and is updated at every time-step. L C0 i is its relaxed length. In a similar manner of the invadopodia membrane, the elastic force at the i-th nuclear membrane (module CN in Table 1),F N E,i , is obtained by using the virtual work theory in structural mechanics. Two types of total elastic energy are considered in the nuclear membrane. Then, F N E,i can be obtained by differentiating the two kinds of total energy, www.nature.com/scientificreports/ where L N i is the length of the i-th line of the nuclear membrane mesh and is updated at every time-step. L N0 i is its relaxed length. κ N L is effective stiffness constants of the line elements of the nuclear membrane (5.0 × 10 -3 N/m) 59,60 . A N i is the i-th mesh area of the nuclear membrane and A N0 i is its relaxed values. κ N A is an effective stiffness constant of area elements of the nuclear membrane (1.0 × 10 -4 N/m 2 ) 61 . F N L,i is the actin crosslinking force by the l-th filamin or α-actinin which is connected to the i-th node in the nuclear membrane, and it is expressed as Computation of viscoelastic cancer cell model. Cancer invadopodia protrusion 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 invadopodia membrane, FTL, ACL, and branched actin network, and double layers of nuclear membrane (PAL and NMS) were numerically coupled to solve for unknown variables associated with node position vectors for all of computational domains (see Table 1). Thereby, all of five submodules were coupled with viscoelastic kelvin-voigt model ( Fig. 2b and 3). To solve five coupled ordinary differential equations in the five domains numerically, these equations should be converted with respect to vectors as followings: where [C] ∈ R 3N cell ×3N cell is dissipation coefficients matrix (Table 1), N cell = N I + N T + N C + N A + N N . Instead of obtaining a solution of Eq. (29) implicitly, we solve Eq. (30) explicitly based two reasons: 1) implicitly solving Eq. (29) is expensive because matrix [C] is required to update at every iteration as the geometric structure of branched actin network is time-varying, and 2) matrix [C] becomes singular without the inclusion of viscous drag coefficients in Table 1. Thereby, explicit forms of five equations in the cell module can be expressed as followings: where F I0 D,i =   Discrete ECM dynamics and reaction-diffusion model. The detailed equations that govern discrete ECM dynamics and reaction-diffusion mass transfer are summarized in S1 Text 21,22 . Visualization of simulated results. All the images (Figs. 1, 2c,d, 4a,b, 5d, 6a,b, 7a,c, S1b, S3 and S4a-c) and supplementary movies (Movies 1, 2, 4, 5, 6, 7 and 8) were created using a commercial software, Tecplot 360 (version 2018 R2, https:// www. tecpl ot. com/ produ cts/ tecpl ot-360/). Cell culture and Preparation of Microfluidic Devices. GFP-expressing MDA-MB-231 breast carcinoma cells were cultured in DMEM. DMEM was supplemented with 10% FBS and 100 U/mL penicillin/streptomycin. Cells were grown in T25 flasks (Thermo Fisher Scientific) until 90% confluent. Then they were removed from the flask surface by incubation with trypsin/EDTA (Gibco) for 2 min, centrifuged, and then resuspended in a 2.5 mg/mL collagen I solution, before seeding into the 3D microfluidic devices.
Preparation of Microfluidic Devices and Invadopodia Imaging. Preparation of 3D microfluidic devices, followed by cell seeding, was performed as previously described in detail 67 . In brief, microfluidic devices were prepared using polydimethylsiloxane (PDMS) and bonded to a glass coverslip. Microfluidic channels were coated with poly-D-lysine (PDL) (Sigma). Before injection of the gel solution into the devices, cells were mixed with the gel solution from the original stock. Following that, the central gel channel, in which invadopodia dynamics were later imaged, was filled with the collagen I solution, at a density of 2.5 mg/mL. Devices were then incubated for 30 min at 37 °C to polymerize. The concentration of NaOH in the collagen I solution was varied to allow accurate control of pH 8.0 during polymerization of the solution. Following the polymerization process, cancer cells were completely surrounded by collagen type I ECM with estimated pore size (1 μm) 68 . After overnight incubation, the microfluidic device was transferred to a confocal laser-scanning microscope (Zeiss) fitted with an environmental chamber operating at 37 °C and 5% CO 2 . Time-lapse microscopy was employed to record cancer cell movement in the 3D collagen I ECM. Protrusive and retractile invadopodia of the GFP-expressing HUVECs (green), and collagen fibers (confocal reflectance microscopy) were simultaneously imaged with a time interval of 3 min for 1.5 h.