Balance of mechanical forces drives endothelial gap formation and may facilitate cancer and immune-cell extravasation

The formation of gaps in the endothelium is a crucial process underlying both cancer and immune cell extravasation, contributing to the functioning of the immune system during infection, the unfavorable development of chronic inflammation and tumor metastasis. Here, we present a stochastic-mechanical multiscale model of an endothelial cell monolayer and show that the dynamic nature of the endothelium leads to spontaneous gap formation, even without intervention from the transmigrating cells. These gaps preferentially appear at the vertices between three endothelial cells, as opposed to the border between two cells. We quantify the frequency and lifetime of these gaps, and validate our predictions experimentally. Interestingly, we find experimentally that cancer cells also preferentially extravasate at vertices, even when they first arrest on borders. This suggests that extravasating cells, rather than initially signaling to the endothelium, might exploit the autonomously forming gaps in the endothelium to initiate transmigration.


Introduction
Immune and cancer cells alike are characterized by their ability to migrate within the vasculature and then to leave the vasculature into different tissues. These processes are crucial for a functioning immune system to fight acute infections [1] or participate in wound healing [2]. However, chronic inflammation or tumor metastases are ultimately also initiated by extravasating immune or cancer cells, respectively [3,4,5]. Hence, while extravasation is critical to cure communicable diseases, it is also a critical contributor to virtually all non-communicable disease, ranging from cancer to asthma, atherosclerosis, rheumatoid arthritis and heart diseases [6,4,7].
Much of the research on extravasation (often termed diapedesis in the context of immune cells) has focused on the role of the extravasating cell during this process, and how it interacts with the endothelial cells of the vasculature through which it is transmigrating. First, the extravasating cell needs to arrest in the vasculature. This may occur through single cells or clusters getting physically stuck in small capillaries, through the formation of adhesions, or both [8,9,10,11,12]. Such adhesion is mediated by molecules including P-and E-selectin, ICAM, VCAM or integrins [13]. The actual process of transmigration can occur through a single endothelial cell (transcellular extravasation) or, more commonly, in between two or more endothelial cells (paracellular extravasation) [14,1].
During paracellular extravasation, it was investigated how the extravasating cell signals to the endothelial cells, leading to weakening of VE-cadherin-mediated cell-cell junctions and subsequently gap formation, through which the cells can transmigrate [8,1]. Gap formation may, for instance, be stimulated by thrombin [15]. As such, molecular signaling events are firmly established as important contributors to extravasation of immune cells. However, on a fundamental level, all the processes involved in extravasation are mechanical processes. Transmigration, like other forms of cell migration, involves the generation of mechanical forces through the actomyosin cytoskeleton [16]. Moreover, the mechanical properties of the endothelium provide passive mechanical resistance [16]. For instance, increased endothelial cell and junctional stiffness will reduce paracellular extravasation rates [17,14]. Interestingly, recent research established that active mechanical properties of the endothelial cells are also critical during endothelial gap formation [18,19,20], and the rearrangements of cytoskeletal structures are associated with changes in barrier function. For instance, a rich actin cortex parallel to cell-cell borders is associated with stabilized VE-cadherin junctions and thus tight barriers [21,22], whereas actomyosin stress fibers pulling radially on junctions can lead to junctional remodeling [18,23]. Additionally, actin-rich pores can actively contract to prevent leakage during extravasation [24]. However, there is still a lack in mechanistic and systems level understanding of the different roles of active and passive mechanical properties of the endothelium.
Mathematical multiscale models are powerful tools to investigate the interplay of different physical drivers of biological processes. Many different approaches have been employed to model and understand the dynamics of epithelial monolayers. Agent based models, where individual cells are explicitly taken into account, include center based models (CBM) [25], vertex models [26,27] and deformable models (DFM) [28,29]. However, these models do not explicitly model cell-cell adhesion dynamics in a way that leads to the experimentally observed gap formation in monolayers of endothelial cells, and can thus not easily be employed to study this problem so crucial for cancer and immune transmigration.
In this paper, we introduce a mathematical multiscale model of the mechanics of an endothelial monolayer where each endothelial cell contains contractile actin structures that may contract radially or in parallel to the plasma membrane. Then, cells are tethered to neighboring cells by cell-cell junctions that can dynamically form and break in a force-dependent manner. We employ this model to investigate the mechanisms of gap formation in an endothelial monolayer. Interestingly, we find that gaps open dynamically in the absence of any extravasating cells. These gaps form preferentially at the vertices where three or more endothelial cells meet, as opposed to the borders in between two cells. This is in line with our experimental data obtained in-vitro from quantifying gap formation of monolayers of human umbilical vein endothelial cells (HUVECs) seeded on glass. We quantify the frequency of gap openings as well as the duration of gap openings, obtaining good agreement between numerical predictions and experiments. Moreover, through multi-dimensional parameter studies, the mathematical model is able to give us insights into the physical and molecular drivers of the gap formation and gap dynamics. The model predicts that active and passive mechanical forces play an important role in the initial gap formation and in controlling size and lifetime of gaps once they initially formed. The catch bond nature of the cell-cell adhesion complexes as well as the forcedependent reinforcement of adhesion clusters may both stabilize junctions in response to forces acting on them. However, while the catch bonds ultimately weaken when forces are increased beyond the maximal lifetime of a single molecular bond, the force-dependent reinforcement will increase adhesion strength with increasing force [30,18]. While the catch bond nature and the force dependence of the adhesion clustering processes both crucially influence gap opening frequencies, we find that gap lifetime and gap size are even more sensitive to the passive mechanical properties of the cell. Increased stiffness of the membrane/cortex and, even more notably, of the actin stress fibers will reduce lifetime and size, since the cells will then increasingly resist opening gaps through counteracting forces. On the other hand, we find that changes in bending stiffness of the membrane/cortex may have gap promoting or inhibiting effects.
Our model predictions of gap opening frequency and lifetime at both cell vertices and borders are validated by experiments observing such gaps in endothelial monolayers in the absence of any extravasating cell. The results thus challenge the paradigm that all extravasating cells primarily cause gap opening through interactions with the endothelium [1,8,31]. We then show experimentally that extravasating cancer cells indeed primarily extravasate at vertices, in line with similar observations for neutrophils [32]. Moreover, we show that cancer cells prefer to extravasate at vertices even when they initially attached to the endothelium at twocell borders. This suggests that, even though extravasating cells can actively interact with the endothelium during transmigration, as shown in earlier studies, they may also take advantage of the autonomous occurrence of a gap, as predicted and verified to occur in our model. In summary, our work highlights the importance of taking the dynamic and autonomous mechanical properties of the endothelium into account when trying to understand gap formation and extravasation. mechanical state is determined by radial contractile actin stress fibers and the cell membrane together with the actin cortex. For simplicity, we combined membrane and cortex into single viscoelastic elements, composed of an elastic spring and a viscous damper, that we refer to, from now on, as membrane elements. The radial stress fibers are also modeled by viscoelastic elements with different mechanical properties from the membrane, similar to a model of epithelial cells [28] (see Fig 1A). Neighboring cells may form cell-cell adhesions at adjacent nodes, and the resulting adhesion bond is modeled through a spring. The passive mechanical properties of the monolayer are thus modeled through a network of connected elastic and viscoelastic elements, similar to models of epithelial sheets [29,28]. Since we are interested in studying the opening dynamics of gaps in the endothelial barrier, we explicitly simulate the dynamical binding of adhesion complexes. Contractions represent myosin motor activity that is known to exhibit randomness [33], so we employ Monte-Carlo simulations to estimate the occurrence of such forces as well as that of protrusive forces due to actin polymerization. The forces are then redistributed across the network of connected viscoelastic elements. Cell-cell adhesion complexes that mechanically link neighboring cells can dynamically bind and unbind in a force-dependent manner. The adhesion complexes in the model provide an effective description of both bonds of cell-cell adhesion molecules (such as VE-cadherin) and bonds of these adhesion molecules to the cytoskeleton. Cadherins and adhesion-cytoskeleton bonds are known to increase their binding strength in response to smaller forces, before they ultimately rupture [34]. This catch-bond type behavior is included in our model, and unbinding is thus simulated through a force-dependent Monte-Carlo simulation. Moreover, the number of VEcadherins in an adhesion complex is modeled through a force-dependent adhesion clustering mechanism, as described in [18,23,35,36,37]. A more detailed description of the mathematical model and its numerical implementation is given in S1 Text.
We employ our endothelial monolayer model to explore the dynamics of endothelial cell junctions. We predict the frequency, size and duration of gaps, as well as the preferred geometrical locations of the gap formation, and compare the predictions with our experimental measurements. The parameters used in the simulations are detailed in S1 Table. After comparing our predictions with the experimental results, we perform sensitivity analyses to investigate how cell mechanical properties, cell-cell adhesion characteristics and myosin generated forces regulate the formation, lifetime and size of gaps in the endothelium.

Summary of major model parameters
Here we present a summary of the major parameters of the model that had a significant impact on our model behavior, and were consequently thoroughly investigated through sensitivity analysis in the remainder of this paper. Table 1 lists all these parameters, and for a complete list and discussion see the Supporting Information. The main parameters investigated are related to cell mechanical properties, adhesion properties or myosin force generated processes.
Cell mechanical properties are dictated by stress fiber stiffness (K sf ), membrane stiffness (K memb ) and bending stiffness (incorporated through a rotational spring constant, K bend ). Stress fiber stiffness controls the rigidity of the interior of the cell, whereas membrane stiffness controls the rigidity of the membrane and the adjacent actin cortex. Bending stiffness acts on the membrane nodes depending on the relative orientation between the edges connecting at a given node.
Adhesion properties are controlled by the mechanical properties of the adhesion complexes and their binding and unbinding rates. Adhesion complex mechanics are modeled by linear springs, controlled by their stiffness constant, K 0 adh . The binding rate depends on distance and can be controlled by the adhesion complex density, ρ adh . We then model the reinforcement of a bond that is already formed by the additional recruitment of adhesive proteins into the bond. Reinforcement is force dependent and can be controlled by the binding rate constant for adhesion reinforcement, k 0 reinf . Unbinding follows a catch bond behavior. The catch bond unbinding curve can be modified through two rate coefficients: k 0 s , which represents a slip bond, and k 0 c , which is the additional parameter characterizing the initial increase in the bond lifetime with force (see S12 Fig).
Then, the model includes contractile forces due to myosin motor activity, and protrusive forces that may arise due to actin polymerization. These forces can be directed radially (following the stress fibers direction) or in a tangential direction (following membrane segments). In the sensitivity analysis we have varied the magnitude of contraction forces in the radial direction (F Radial ) and in the tangential direction (F cortex ).

Gaps open preferentially at vertices
Fig 1B and 1C and S1 Movie show typical simulations of the monolayer dynamics of the computational model. We observe that gaps open preferentially at vertices, i.e. the intersections of three or more cells, as opposed to the border between two cells. We have quantified this by counting the total number of gaps formed as well as their lifetime at borders and vertices of the cell in the center of the monolayer, and showed that our model predictions are in line with the experimental observations ( Fig 1G and 1H). These experiments were performed by seeding HUVEC cells on glass, where they formed a continuous monolayer. The gaps were experimentally quantified through inspection of visible gaps within the VE-cadherin-GFP signal in the monolayer (arrows in Fig 1E and 1F). Controls simultaneously showing VE-cadherin-GFP and CD31 staining show that the VE-cadherin gaps are also visible in the CD31 staining, indicating that the VE-cadherin gaps correspond to real physical gaps between two or more cells S7 Fig (see Methods for further details of the experimental setup and quantification). Vertices are points where more than two cells exert forces and where tangential force components naturally propagate to. Therefore, it is expected that stress concentrates at the three cell vertex rather than at the two cells borders, and the simulations confirm this hypothesis ( Supplementary S9 Fig and S3 Movie). The forces on adhesion clusters at the vertices are thus more likely to exceed the corresponding force of maximal lifetime of the bonds, as will be discussed in more detail below.

Mechanical properties of cell-cell adhesion complexes limit endothelial gap opening frequency
We study how variations in the mechanical properties of the cells, the cell-cell adhesion complexes or force variations affect the rate of gap formation. Fig 2A and 2B show how passive mechanical properties of the cell affect both the frequency (Fig 2A) and the location of the gap openings ( Fig 2B). Increasing stiffness of either the membrane or the stress fibers provokes a decrement of the gap generation frequency (Fig 2A and S4 and S5 Movies). This is intuitive, since increasing stiffness stabilizes the movements of cells and makes the monolayer less dynamic. On the other hand, the location of the gap openings (i.e. whether they occur at a vertex or border) is critically affected by membrane stiffness at low values, until it stabilizes for intermediate and high membrane stiffness. In contrast, stress fiber stiffness affects gap location for very high stiffness, where gaps are almost fully prevented from opening at the borders ( Fig  2B). Interestingly, increasing bending stiffness first increases gap generation up to a maximum point, before it leads to a decrease in gap opening frequency (Fig 2A). For small to intermediate bending stiffness, the frequency of gap openings increases, since bending stiffness is critical for effective force propagation between neighboring adhesion sites at vertices. When a single adhesion complex ruptures, bending stiffness leads to increased forces on neighboring adhesion complexes. After a peak in gap opening frequency at intermediate bending stiffness, a drop in the gap formation is observed for higher bending stiffness. This is caused by the resulting stabilization of the existing gaps at vertices. This high bending stiffness opposes sharp corners of the membrane at vertices and thus favors stable gaps that are permanently open, implying no new gaps are formed (S7 Movie). On the other hand, at cell borders, a high bending stiffness implies that if a single adhesion cluster is ruptured, the forces on it are redistributed across many neighboring adhesion sites and this stabilizes the borders ( Fig 2B).
Turning to the role of cell-cell adhesion complex properties, our model shows that as the junctions become more stable, gaps open less frequently ( Fig 2D). To increase cell-cell junction stability, we increase the mechanical stiffness of individual adhesion bonds, or the density of adhesion molecules. These results are in line with previous experimental work [14], which reported that more stable cell-cell junctions result in fewer transmigrating cells. While the total number of gaps at either vertex or border decreases with increasing cell-cell adhesion complex stiffness or cell-cell adhesion density available for binding, we see that there are no significant differences between gaps generated at the vertex and gaps generated at the borders ( Fig 2E). Fig 2G and 2H show the impact of changing the cortical and radial forces, where the total force is kept constant (when the radial force decreases, the cortical force is increased by the same magnitude). This is biologically relevant since cells are known to shift their cytoskeletal compartments in a context dependent manner [38]. In fact, cell monolayers subjected to shear flow have been reported to increase cortical actin while decreasing stress fibers [14]. The second column (B, E, H) shows the ratio of gaps that occur at a two cell border to the gaps that originate at a three cell vertex. Error bars show to the standard error. The third column (C, F, I) shows the impact of a two parameter variation in the gap opening frequency. The first row shows results from varying cell mechanical properties (stress fiber, membrane and bending stiffness (A, B)), and simultaneous variations of stress fiber and membrane stiffness (C). In the second row (C-F) properties of cell-cell junction are changed: adhesion stiffness and adhesion density. The third row shows results for increasing cortical force, keeping the total force constant (G and H) or varying both forces simultaneously (I). https://doi.org/10.1371/journal.pcbi.1006395.g002 Endothelial cells in particular, are known to exhibit both radial and tangential stress fibers with a different effect on gap opening dynamics [39]. As the force shifts from radial to cortical forces, total gap formation fluctuates with a slight increase as cortical forces increase ( Fig 2G). For high cortical forces, the gaps also clearly tend to localize more at the vertices ( Fig 2H). This is because contractions parallel to the membrane result in force concentrations at the vertices. For very high cortex forces, the typical stresses on adhesion clusters at the vertices may thus be higher than the force where the lifetime of catch bonds peaks (Supplementary S12 Fig), explaining the small increase in the number of gaps formed (Fig 2G). On the other hand, we will later show that these gaps formed at high cortical forces are typically small and have a short lifetime, limiting their potential for extravasation (see Fig 3I and 3J).
To take into account that molecular or physical perturbations may simultaneously affect multiple parameters, we now study how variations of pairs of these parameters at the same time may influence the monolayer integrity and the localization of the gap formation. Although, we have previously seen in Fig 2A that membrane and stress fiber stiffness have a similar effect on the gap opening frequency, in Fig 2C we can observe how the effect of varying stress fiber stiffness is clearly predominant over the effect of varying membrane stiffness. Fig  2F shows the impact of varying cell-cell adhesion stiffness and cell-cell adhesion complex density available for binding. Interestingly, there is a synergy between both parameters on regulating gap opening frequency, as evident through the curved shape of the levels of equal gap opening frequency (Fig 2F). In Fig 2I we show the combined role of cortex and radial forces, thus not keeping total force fixed as in Fig 2G and 2H. This confirms that total force is the main driver of gap opening frequency, as opposed to a redistribution of forces between cortex and stress fibers (Fig 2I).

Passive cell-mechanical properties dictate endothelial gap lifetime and size
The lifetime and size of a gap are physical parameters that may also limit a cancer or immune cell's potential to extravasate through the monolayer. Here, we show how the lifetime and size of a gap are influenced by cell mechanical and junction properties, without the presence of extravasating cells (Fig 3). We observe that membrane stiffness has a marginal influence on the life time of the gap, whereas increasing stress fiber stiffness clearly reduces the time that a gap is open and the gap size (Fig 3A and 3B). Indeed, higher stress fiber stiffness will result in mechanical resistance to an opening gap and thus inhibit the propagation of the defect in the cell-cell junctions, leading to a stabilization of the monolayer (see S4 and S5 Movies). The dominance of stress fiber stiffness over membrane stiffness in regulating lifetime and size remains valid in a broad range of parameter values (Fig 3C and 3D).
Interestingly, increasing bending stiffness to high values may increase gap lifetime (Fig 3A). This is because higher bending stiffness will resist deviations from straight membranes. Thus, at straight borders, higher bending stiffness will resist gap openings whereas at vertices with high curvature, cells are more likely to adapt their shape resisting high curvature, thus favoring opened gaps. The dynamics of the monolayer for low bending stiffness is shown in S6 Movie. Fig 3E and 3G show that adhesion complex stiffness and density at low values do not have a big impact on lifetime, however as they increase, lifetime starts to decrease. Both stiffness and density have a similar effect, since the total stiffness of an adhesion complex depends on both density and single bond stiffness (Eq. S10). Higher stiffness of the adhesion complex leads to more passive mechanical resistance to gap openings, and this effect dominates for high stiffness. The level of noise due to repeats of our MC simulations is higher for these adhesion parameters than for the parameters determining cell mechanics. Likewise, for the gap size, the stabilizing effect of both adhesion complex stiffness and density dominates and leads to a reduction in gap size (Fig 3F and 3H). However, the effect of increasing the density is slightly stronger than that of increasing single bond stiffness. This is because the density affects not only adhesion complex stiffness (Eq. S10), but also the rate of forming new adhesion complexes (Eq. S9) and the rate of reinforcing existing bonds (Eq. S11). These effects together thus synergize to stabilize gaps and prevent them from growing too large.
Earlier, we have shown that a shift in the force (from radial to cortical) produces an increment in gap formation (Fig 2G). Fig 3I and 3J show that this shift in the force reduces gap lifetime and size. This indicates that, although the frequency of opening is increased, these gaps are smaller and last shorter in time which may reduce paracellular extravasation, as suggested in previous experimental work [14]. Combined changes of cortical and radial force show that although both kinds of forces are needed to increases gap size and lifetime, the impact of radial forces is clearly predominant over the impact of cortex forces (Fig 3K and 3L). This is intuitive, since radial forces clearly separate cell borders generating bigger gaps and make them harder to close, whereas cortical forces distribute forces to vertex regions. This does not provoke large cell deformations, which is reflected in the low impact on the gap size and lifetime observed.

Catch bonds facilitate regimes of maximal endothelial stability
In Fig 4A-4C we show the impact of varying the catch-bond unbinding parameter k 0 c that shifts the location of the peak of maximal lifetime of a single catch bond, while we maintain the actual maximum value through simultaneously shifting the slip-bond unbinding parameter k 0 s (Eq. S12 and S12 Fig). We observe that for a pure slip bond (corresponding to k 0 c ¼ 0), gaps occur at a higher rates than for small nonzero values of k 0 c . Increasing k 0 c further leads to a minimum in gap opening frequency, from which the frequency increases again. This minimum corresponds to a maximum of stability, where forces on the adhesion complexes are similar in magnitude to the peak of stability of the catch bond. Consequently, shifting the location of that peak even further towards higher forces (by increasing k 0 c even further) means we destabilize the catch bonds again. Note that the gap lifetime and size of gaps are much less influenced by the location of the catch bond maximum than the gap opening frequency.
In Supplementary S11 Fig, we show histograms of the forces on adhesions comparing the number of bound clutches, the number of unbinding events, and the ratio of unbound to total bonds for slip bonds (k 0 c ¼ 0) to the catch bond with reference values (k 0 c ¼ 0:27s À 1 ). We see that adhesions in the catch bond case bear and disengage at higher forces than for the slip bond case, confirming that the typical forces on bonds are of such magnitude that the catch bond nature stabilizes the junctions.
In Fig 4D-4F we modify the reinforcement binding rate k 0 reinf to check the influence of the reinforcement. This is different from the previous analysis where the adhesion complex density available for binding was changed, since now the binding probability based on distance is not affected (Eq. S9). However, we see the same trend of increasing stability with increasing k 0 reinf (Fig 4D), in line with the result obtained from varying cadherin density (Fig 2D), suggesting that binding is mainly regulated by this reinforcement process. Similar to the catch bond, we see that adhesion reinforcement is less important in determining gap size or lifetime ( Fig  4E and 4F) than in regulating gap opening frequency.

Force fluctuations and distribution regulate gap opening dynamics
We have shown that both the magnitude of forces and the cytoskeletal compartment that generates the forces (stress fibers or cortex) affect gap opening frequency, size and and lifetime. Besides these broad compartments, many other biological and physical parameters affect how forces ultimately act on cell-cell junctions: Forces may act in a directed manner due to larger parallel actin bundles and synchronous myosin activation, e.g. initiated through waves of activators [15], or may act more randomly [33]. We test variations in force applications through parameters that affect the transition time when forces change (t Force Transition ), through spatial force distributions and through the velocity at which forces are modified. In Fig 4G we

observe how increasing the force transition time t Force
Transition slowly reduces the gap opening. This is due to the fact that a slower, persistent application of forces leads to a redistribution of the forces through rearrangement and remodeling of the cell. It is consistent with experimental works that showed that force fluctuations influence gap opening dynamics [15].
Then, distributing the same radial forces over several adjacent stress fibers reduces gap opening frequency (Fig 4H). More spatially distributed forces are less capable of damaging cell-cell junctions than localized peak forces, since such high peak forces are required to overcome the catch bond maximal lifetime. Likewise, high peak forces lead to longer lifetime and larger size of the resulting gaps (Supplementary S13C and S13D Fig).
Next we observe the effect of force persistence in time. We vary the force recalculation time parameter (equally for all forces) in Fig 4I. Results show that the time that forces are applied does not have a big influence on gap formation. This suggests that cells are able to adapt to forces in longer time scales and therefore it is not the time that forces are applied what regulates gap formation, but the transitions of force fluctuations and their spatial distribution.

Cancer extravasation mimics autonomously occurring endothelial gap dynamics
To demonstrate that the geometry of the gap opening dynamics is physiologically relevant, we quantified the characteristics of extravasating cancer cells through monolayers of HUVECs, as shown in S9 Movie. Here, a tumor cell is seen transmigrating through an endothelial monolayer at a tricellular junction as delineated by VE-cadherin GFP, followed by gap-closure after the tumor cell has completely cleared the barrier. Fig 5A shows the ratio of tumor cells that extravasated at vertices, relative to borders. We see that tumor cells preferentially extravasate at the vertices, in line with the previously observed increased frequency of gaps opening there ( Fig 1G) and similar observations of extravasating neutrophils [32]. Moreover, even if cancer cells initially arrest at the border between two endothelial cells, they are much more likely to extravasate at a vertex at later points in time, rather than at the border where they initially attached to, perhaps first through migration on the surface of the endothelium and subsequent preferential attachment to points of exposed basement membrane as a result of inherent EC junctional dynamics (Fig 5B). This could suggest that in addition to the possibility of cancer cells actively signaling to open gaps in the endothelium, endothelial barrier dynamics itself can also present the cancer cells with opportunities to begin the transmigration process.

Discussion
The computational model presented in this paper allowed us to study how gaps in an endothelial monolayer initially open, grow, stabilize and finally close, and we identified which physical properties dominantly regulate each stage.
The model simulates a cell monolayer in two dimensions. Adhesion between cells is simulated through binding or unbinding of adhesion complexes located on adjacent cells. These adhesion complexes are dynamically engaging and disengaging as the myosin generated forces cause cell deformations. Because of cell-cell adhesion rupture, gaps between the cells are formed. To perform our simulations, the model is based on a number of assumptions that simplified the model. First of all, due to the typically small height (about 3μm) of endothelial cells [40], we neglected the third dimension perpendicular to the monolayer. However, disruptions of the spatio-temporal dynamics of adhesion molecules and cytoskeletal organization in the third dimensions are likely to impact gap formation. Incorporating such effects into our model would consequently require a 3D model of a cell with more detailed descriptions of the subcellular mechanics. However, the purpose of our model was to demonstrate the broad impact of subcellular mechanical structures on gap formation. For this reason, we modeled the cells in two dimensions and included only radial stress fibers and contractile actin fibers parallel to the membrane. This was motivated by experiments that indicated different roles of these actin structures on gap formation [23]. Each discrete adhesion complex is simulated as one cluster to simulate recruitment of proteins such as vinculin or talin, without increasing the total number of components in the simulation. Myosin generated forces included in the model are assumed to occur only in the direction of the stress fibers or membrane.
Supplementary S14 Fig summarizes some of our key conclusions: By comparing our reference case with extreme variations of very low stress fiber or membrane stiffness, we see that both passive mechanical properties and adhesion complex properties are important in controlling gap opening frequency (Supplementary S14A Fig). On the other hand, the lifetime and especially size of gaps increases significantly with lower stress fiber or membrane stiffness, since the softer cells are more likely to deform and adapt in response to the opened gap (Supplementary S14B and S14C Fig). We also verify that stress fiber stiffness influence is stronger than membrane stiffness influence. In contrast, properties of the cell-cell adhesions strongly affect the frequency of the gap openings, but less so their lifetime or size. Indeed, decreasing the density of adhesion bonds or the adhesion stiffness strongly increase the frequency of forming gaps (Supplementary S14A Fig), while only marginally affecting the size or lifetime of the gaps (Supplementary S14B and S14C Fig). These data thus summarizes our biological model where adhesion properties control the initial formation of gaps, while cell mechanical properties are critical in limiting the size and duration of opened gaps.
Our results that gaps open more frequently at vertices than at borders were true over wide ranges of parameters (Fig 2B, 2E and 2H). Only extremely small bending stiffness led to similar frequencies of gaps at vertices and at borders (Fig 2B). These results also show that earlier experimental observations, where neutrophils were found to extravasate preferentially at endothelial cell vertices, [32], can be explained through the mechanical dynamics of the endothelial monolayer alone. Consequently, this may be a general mechanism for extravasating cells, and we found a similar behavior with extravasating cancer cells (Fig 5). This finding is complementary to the extensive literature that suggests that chemical or mechanical signaling of extravasating immune or cancer cells to the endothelium facilitates extravasation [1,8,31]. There are many potential hypotheses why both the autonomous dynamics of the endothelial monolayer and the bidirectional signaling with the extravasating cells may play a role during extravasation: It may be that initial autonomously forming gaps are important for extravasating cells to sense a gap and they consequently signal to widen the gap or to keep it open. The gap sizes that the model predicts are of the order of magnitude of a few microns, which is enough for extravasating cells to protrude through the gap. Our previous study indicates that tumor cells can squeeze significantly when transmigrating through artificial gaps [16], so the autonomous gaps may be of sufficient size for complete transmigration. Nevertheless, endothelial gaps may widen during transmigration, so crosstalk between the transmigrating cell and the endothelium likely remains an important factor contributing to the likelihood and speed of extravasation. Then, whether bidirectional signaling or autonomous gap formation dominates the process may be cell type specific. For instance, it is still a major research question why certain cancer cells preferentially metastasize to certain organs [41]. We may speculate that not only the signaling of the specific primary tumor cells with an organ-specific type of endothelial cells influences the likelihood of extravasation [8]. Also, the mechanical properties of the endothelium of the target organs will likely play a major role. Our flexible modeling framework was tested with a HUVECs monolayer, yet, by changing the physical parameters of the model, it may be quickly adapted to other endothelia. Besides testing our model with different endothelial cells, some other important steps towards validating our model conclusions in vivo will be to test our model with more realistic three dimensional microvasculature with blood flow, embedded in extracellular matrix and surrounded by supporting cells such as pericytes, fibroblasts or, for brain, astrocytes [41,1]. Such real, in vivo microvasculature consists of vessels that are curved and exposed to shear stresses due to the flow. That, in turn, may be affected by extravasating cells that may obstruct blood flow. Similarly, matrix stiffness was shown to affect endothelial monolayer integrity [42,43]. Some complications in validating our results in vivo involve the lack of available in vitro cultures that are required to provide high throughput, microscopy resolution and level of experimental control that is lacking in vivo, making direct comparison of computational models to in vivo experiments unfeasible. However, the recent rapid progress in developing more complex and organ specific in vitro assays of 3D microvasculature will make such validations feasible in the near future [44,45,46].
Our model is based on a number of simplifications. We do not consider the effect of extracellular matrix and substrate stiffness properties on monolayer integrity, despite the known effect of these properties on cell mechanics. It is important to remark that cells on glass may behave very differently than in vivo endothelial vessels. Our model also does not include the effect of fluid pressure or tangential stress due to fluid flow. Pressure and blood flow would induce additional forces over the monolayer that could affect gap generation processes. For example, it was observed [14] that tangential flow could induce the strengthening of cell-cell junctions, therefore reducing paracellular extravasation.
Modeling such complex environments presents a great challenge to both in vitro and in silico models. It is therefore essential to justify assumptions that can reduce this complexity and make the model development feasible. Here, we have assumed that the mechanics of the inside of a cell is determined by a fixed number of stress fibers, although it is known that inside the cell there are different polymer structures such as microtubules and intermediate filaments. Moreover, actin filaments are not fixed in time but appear and disappear depending on their stability and polymerization rates. To simulate all of this with high accuracy would require a completely different model in which the computational cost that would exceed current capabilities. For the purpose of this project, we focused on incorporating essential cell mechanical structures that have been implicated in the regulation of gap formation, and modeled a fixed number of stress fiber similar to other works [28,29]. Similarly, we have simulated adhesion complexes as discrete elements that can bind two membrane points of neighboring cells. In real cells, adhesion complexes between cells are formed by a great variety of proteins such as VE-cadherins, α-catenin, talin or vinculin. While the spatio-temporal dynamics of each of these adhesion molecules likely influences gap formation, no computational model can currently explain their precise organization in adhesion complexes and their resulting effect on gap formation. Consequently, our model included an effective term that describes the force dependent recruitment of adhesions, as observed in different experimental studies [23,20].
Moreover, there are also challenges to the mathematical modeling of complex 3D microvasculature. Modeling of epithelial sheets in 3D has proved challenging, with some recent interesting progress after decades of mainly focusing on epithelial monolayers in 2D [47,48,49]. These models are based on frameworks such as vertex models, where the dynamics of each cell is incorporated into the dynamics of vertices between cells. There are many other modeling frameworks that can capture different aspects of the complex cell behavior, such as cell based models [50], immersed boundary models [51] or subcellular element models [52,53]. These modeling frameworks are, however, not directly suitable to predict the formation of gaps at either vertices or borders. Given these challenges, is was paramount to establish a 2D mathematical model of an endothelial monolayer that was validated with novel experiments and that was able to lead to insights into the mechanisms of endothelial gap formation.

Generation of HUVEC monolayers and image acquisition
Human umbilical chord vein cells (HUVECs) were transduced with VE-cadherin-GFP using methods described previously [45]. HUVECs at P7-10 were seeded onto 35 mm glass bottom Mattek dishes (at 3 × 10 5 cells/dish), which had been plasma treated for 30 seconds previously. Cells were allowed to grow to confluence (beyond 100%) in EGM-2MV (Lonza) for 3 days before imaging. Dishes were imaged on an Olympus FV1000 confocal microscope with magnifications of 30X (oil immersion), under an environmental chamber set at 37C and 5% CO2. The chamber was equilibrated for * 30 min prior to the start of image acquisition. For timelapse videos of junctional dynamics, z-stacks of 40μm (4μm steps) were taken at intervals of 3 minutes.

Analysis of junctional disruption dynamics
Time-lapse images were appended and analyzed manually on ImageJ. A single unique junctional disruption is defined as a vertex or border with an observed gap of greater or equal than 3μm, and are preceded and proceeded at some point in time with a closure (no visible gap in fluorescence greater than 0.6μm). The number of junctional disruption events was counted for each border and vertex of an image over a total time period of 2 hours. Vertices and borders belonging to the same cells were still considered to be unique.

Analysis of tumor cell extravasation
Tumor cells were suspended in EGM-2MV (Lonza) and a concentration of 15,000 cells/mL, and 1mL of the suspension was gently added to each HUVEC monolayer. Cells were allowed to settle first for *10 minutes before acquisition of t = 0 images. For quantification of extravasation, z-stacks were taken at 3μm steps at an endpoint of 6 hours to image the entirety of the tumor cell and endothelial monolayer. Any tumor cell that has breached the endothelial layer as evidenced by protrusion extension across and beneath the endothelial layer was considered as "extravasated". Delineation of the endothelial barrier is visualized via CD31 staining (Biolegend, Cat # 303103) for 30 min in EGM-2MV at 37C and 5% CO2 prior to imaging. Gap opening frequency, average lifetime and size in each column. Note that the point where x and y coordinates are 1 corresponds to the reference case. Error bars represent standard error. In the first row (A, B, C), results for medium and dashpot viscosities of the stress fibers and membrane and varied. Increasing viscosity reduces node movement, stabilizing monolayer dynamics. Medium viscosity has a higher effect on gap opening dynamics since it affects the overall timescale of all mechanical parts of the model. The stress fiber dashpot strongly influences gap lifetime and size; this is similar to the dominating effect of stress fiber stiffness over membrane stiffness on gap lifetime and size (Fig 3A and 3B). The second row (D, E, F) shows the effect of varying the constant for remodeling rate. Increasing the remodeling rate implies that cells are able to adapt their permanent shapes faster in response to deformations. Therefore, the frequency of gap openings increases with the remodeling rate (D). The gap lifetime and size broadly also increase, but less strongly then the opening frequency. (TIF) S1 Movie. Simulation of the endothelial monolayer dynamics. Gaps are more likely to appear in the vertex of three cells than at a two cell border. Green denotes the cell membrane, red the inside of a cell, with darker red being the stress fibers. Parameters are the reference values as in S1 Table. (MP4)