Multiscale model of integrin adhesion assembly

The ability of adherent cells to form adhesions is critical to numerous phases of their physiology. The assembly of adhesions is mediated by several types of integrins. These integrins differ in physical properties, including rate of diffusion on the plasma membrane, rapidity of changing conformation from bent to extended, affinity for extracellular matrix ligands, and lifetimes of their ligand-bound states. However, the way in which nanoscale physical properties of integrins ensure proper adhesion assembly remains elusive. We observe experimentally that both β-1 and β-3 integrins localize in nascent adhesions at the cell leading edge. In order to understand how different nanoscale parameters of β-1 and β-3 integrins mediate proper adhesion assembly, we therefore develop a coarse-grained computational model. Results from the model demonstrate that morphology and distribution of nascent adhesions depend on ligand binding affinity and strength of pairwise interactions. Organization of nascent adhesions depends on the relative amounts of integrins with different bond kinetics. Moreover, the model shows that the architecture of an actin filament network does not perturb the total amount of integrin clustering and ligand binding; however, only bundled actin architectures favor adhesion stability and ultimately maturation. Together, our results support the view that cells can finely tune the expression of different integrin types to determine both structural and dynamic properties of adhesions.


Author summary
Integrin-mediated cell adhesions to the extracellular environment contribute to various cell activities and provide cells with vital environmental cues. Cell adhesions are complex structures that emerge from a number of molecular and macromolecular interactions between integrins and cytoplasmic proteins, between integrins and extracellular ligands, and between integrins themselves. How the combination of these interactions regulate adhesions formation remains poorly understood because of limitations in experimental approaches and numerical methods. Here, we develop a multiscale model of adhesion assembly that treats individual integrins and elements from both the cytoplasm and the

Introduction
As the linker between cytoskeletal adhesion proteins and extracellular matrix ligands, integrins play a vital role in the formation of adhesions and profoundly influence different phases of cell physiology, such as spreading, differentiation, changes in shape, migration and stiffness sensing [1][2][3][4][5].
Integrins are large heterodimeric receptors, with a globular headpiece projecting more than 20 nm from the cell membrane, two transmembrane helices, and two short cytoplasmic tails that bind cytoskeleton adhesion proteins (see Fig 1A). In order to form adhesions, integrins undergo lateral diffusion on the cell membrane, switch conformation from bent to extended, and change chemical affinity for extracellular matrix (ECM) ligands, E IL . Integrins also assemble laterally, owing to interactions with talin [6,7], kindlin [8], or glycocalyx [9], and can grow nascent adhesions into mature adhesions [10][11][12]. Integrin diffusion, activation, ligand binding, and clustering occur at the individual protein scale, but their effects can also be reflected on the cellular scale, resulting in a multiscale biological process. Simulations of adhesion assembly based on allatom approaches are too detailed and computationally demanding to capture adhesion formation from multiple integrins. Instead, highly coarse-grained, CG, approaches based on Brownian Dynamics can condense the description of individual proteins into a few interacting CG "beads" that can recapitulate the emergent dynamics of complex biological systems from its individual components (see, e.g., Refs [13][14][15][16][17] for the example of cytoskeleton networks).
Single-protein tracking experiments combined with super-resolution microscopy and computational methods have helped extract physical properties of different integrin types. β-1 and β-3 integrins were found to have diffusion coefficients of 0.1 and 0.3 μm 2 /s, respectively [43]. β-1 integrins also maintain their active conformation longer than β-3 integrins. Freeenergy energy differences between active and inactive states revealed activation rates for β-3 integrins about 10-fold higher that β-1 integrins [44,45]. The intrinsic ligand binding affinity, E IL , of β-1 integrins for soluble fibronectin is about 10-50 fold higher than β-3 integrins, spanning an overall range for the two integrins of 3-9 k B T [46]. β-1 integrins display a catch bond and adhesion strength-reinforcing behavior and are stationary within adhesions [47][48][49]. β-3 integrins, on the other hand, rapidly transit from closed to open conformations, break their bonds from ligands more easily under modicum tensions, and undergo rearward movements within adhesions [26,38]. How these differences in diffusion, rate of activation, ligand binding affinity, and bond dynamics reflect on the assembly of nascent adhesions and on the probability of adhesion maturation remains elusive.
In this paper, we show that mixed populations of β-1 and β-3 integrins localize to both nascent and mature adhesions, suggesting that there could be important interactions between the two types of integrins. To address this question, we have developed a highly CG model of adhesion formation, based on Brownian Dynamics (Fig 2) and study how nanoscale physical properties of different types of integrins interplay in the assembly of nascent adhesions. The CG model treats individual integrins as point particles within an implicit cell membrane and includes actin filaments as explicit semiflexible polymers (Fig 2A). By incorporating nanoscale physical properties of individual integrins, sequential interactions and feedback mechanisms between integrin, ligands and actin filaments (Fig 2B-2D), the model is used to characterize the formation of micrometer-size adhesions at the cell periphery in a multiscale fashion. Our calculations show that integrins with high E IL and enhanced bond lifetimes, such as β-1 integrins, facilitate ligand binding, transmission of traction stress, and engagement of actin networks. By contrast, integrins with low E IL and lower ligand bond lifetimes, such as β-3 integrins, are correlated with clustering, repeated cycles of diffusion and immobilization, and weak engagement of actin filaments. The architecture of actin filaments does not impact the amount of ligand binding and integrin clustering, but determines the probability of adhesions maturation, consistent with previous experimental findings [50]. Collectively, our data reveal important insights into adhesions assembly that are currently very challenging to obtain experimentally. The data supports the general view that cells, by controlling physical nanoscale properties of integrins via expression of specific types, can regulate structural, dynamical, and mechanical properties of adhesions. representation of active, fully extended α IIB β 3 integrin, which is closely related to α v β 3 [82]. (B) Representative images of human foreskin fibroblasts (HFF) fixed after 60 minutes of spreading on fibronectin coated glass coverslips and immunostained for actin and either β-1 or β-3 integrins.

Both β-1 and β-3 integrins localize in nascent and mature adhesions
Motivated by our recent work on integrin catch-bonds regulating cellular stiffness sensing [5], we sought to investigate how interactions between different integrins could affect adhesion formation. Immunostaining in Human Foreskin Fibroblasts (HFF) for actin and either β-1 or β-3 integrins revealed that both types of integrins localize in nascent adhesions at the cell leading edge and in mature adhesions at the end of actin stress fibers ( Fig 1B). This suggests that potential interactions between the different adhesion populations could be important during adhesion formation. To address this question, we developed a computational model to investigate how the nanoscale properties of different integrins affect adhesion formation and stability.

Integrin clustering decreases with ligand binding affinity
Since β-1 and β-3 integrins differ in ligand binding affinity, E IL , and strength of pairwise interactions, E II [44,45], we use the CG model to test how variations in E II and E IL impact adhesion assembly in terms of the amount of integrin clustering, ligand binding, and spatial arrangement of adhesions. Different morphological arrangements of integrin adhesions are detected (Fig 3A-3C). For high E II and low E IL , clustering is promoted (Fig 3D), but only a few integrins are bound to ligands (Fig 3E), resulting in few large integrin clusters ( Fig 3A). Conversely, for low E II and high E IL , only a few integrins cluster ( Fig 3D) while ligand-binding is promoted (Fig 3E), resulting in many ligand-bound integrins and few small integrin clusters ( Fig 3B). When E II and E IL have intermediate values, a mix of big clusters of integrins that are weakly bound to the substrate and smaller, ligand-bound clusters co-exist ( Fig  3C). By systematically varying E II and E IL , morphological regions differing in size and number of ligand-bound integrins versus clusters were precisely identified. A region of few large clusters exists for E II >3 k B T and E IL < 3 k B T; a region of many small adhesions exists for E II < 3 k B T and E IL > 3 k B T; the rest of the parameter space shows co-existence of intermediate-size clusters and ligand-bound integrins (Fig 3D-3E). The fraction of ligandbound integrins increases with E IL and is independent from E II . (Fig 3E). By contrast, clustering is not independent from E IL and is promoted when E IL is low (Fig 3D). In the model, when active, integrins can bind free ligands and cluster, when in close proximity of a ligand or another active integrin, respectively. Since the number of ligands is higher than the number of integrins, the probability for an integrin to find a free ligand is higher than that of finding an active integrin. Therefore, clustering increases less with E II when E IL is high than when E IL is low (Fig 3D). This indicates that integrin clustering and ligand binding are competing mechanisms.
Together, our results show that different arrangements of nascent adhesions can be achieved depending on E II and E IL . When we use high E II and low E IL , as for β-3 integrins, clustering is enhanced, and ligand binding reduced; when we use high E IL and low E II , as for β-1 integrins, clustering is reduced, and ligand-binding promoted. Thus, the competition between clustering and ligand binding can be determined by the integrin type. However, β-1 and β-3 integrins also differ in their rates of activation, which can lead to differences in this competition, by promoting clustering at high E IL . Therefore, we next aimed to understand how activation rates, combined with variations in E II and E IL , impact clustering and ligand binding.

Rate of integrin activation increases clustering and ligand binding
Competition between integrin clustering and ligand binding can be determined by the difference in activation rate between β-1 and β-3 integrins. By varying k a from 0.005 s -1 to 0.5 s -1 , our model shows that both clustering and ligand binding are promoted (Fig 4A and 4B). Using E II = 5 k B T and varying E IL from 3 k B T to 11 k B T, clustering is independent from E IL (Fig 4A), while overall ligand binding increases with E IL (Fig 4B). Clustering is mostly set by the strength of pairwise interactions between integrins, E II . It can be promoted by low E IL. and high k a , leading to a higher number of integrins able to diffuse and cluster ( Fig 4A). Ligand binding is proportional to E IL at all k a . In experiments, variations in integrin activation rate are tied to variations in ligand binding affinity, making it unclear whether it is k a or E IL that determines organization of nascent adhesions. Our model shows that the rate of integrin activation set the level of the competition between ligand binding affinity and strength of pairwise interactions ( Fig 4A).
Experimentally, Mn 2+ or antibodies are typically used to modulate ligand binding affinity [51][52][53][54]. Both of these approaches, however, not only increase ligand binding affinity, but also the lifetime of the ligand bond. The increase of the ligand bond lifetime can be formally represented using a catch-bonds [55], where ligand unbinding rates decrease under tension and promotes stress transmission from the adhesions [56]. Therefore, we next used the model to test how variations in catch bond kinetics, combined with differences in the relative amount of β-1 and β-3 integrins, modulate ligand binding and stress transmission.

Distribution of tension on integrins depends on bond dynamics
Since nascent adhesions transmit tension between the cytoskeleton and the ECM, we next asked how mixing integrins with different load-dependent bond kinetics impacts ligand binding and transmitted tension. The β-1 and β-3 integrins both behave as catch bonds that differ for unloaded and maximum lifetimes ( Fig 2C). In the model, an increase in the percentage of β-1 integrins while keeping the rest as β-3 integrins, increases ligand binding from about 5% to 35% when using actin flow speeds below 15 nm/s (Fig 5A). The percentage of ligand-bound integrins is in direct proportion to the amount of β-1 integrins (Fig 5A). At actin flow speeds below 15 nm/s, traction stress and flow rate are positively correlated, while at higher flows they are inversely correlated (Fig 5B), in agreement with previous findings [60,61]. Interestingly, variations in the relative fractions of the two integrin types do not affect the average tension on each integrin-ligand bond ( Fig 5B). Below 10 nm/s actin flow, the minimum separation between ligand-bound integrins decreases from about 120 to 10 nm by increasing the fraction of β-1 integrins (Fig 4C). Stable adhesions, with minimum separation between ligand-bound integrins of 70nm, form with at least 20% β-1 integrins ( Fig 5C). Together, our results show that the relative fractions of β-1 and β-3 integrins cooperate with actin flow to determine ligand binding and adhesion stability.

Actin architecture does not impact integrin clustering and ligand binding, but changes the physical organization of integrins in adhesions
Interactions of adhesions with a cytoskeleton network play important role in several cell activities, including spreading and migration. The actin cytoskeleton exists in different architectures, depending on the cell location and function. Therefore, we next considered how the architecture of the actin cytoskeleton can impact the formation of adhesions. We incorporated in the model explicit actin filaments, using random, crisscrossed, and bundled architectures (Fig 6A-6C). The model assumes that ligand-bound integrins can interact with actin filaments, and that binding to actin increases integrin activation rate, as detected experimentally [63,64]. Increasing the fraction of β-1 integrins, ligand binding increases independent of network architecture ( Fig 6D). By contrast, integrin clustering remains at about 20-30% when a percentage of β-3 integrins is used. When only β-1 integrins are used, integrin clustering decreases of about 3-fold, independent from network architecture ( Fig 6E). The number of ligandbound integrins with a separation less than 70 nm is enhanced using a bundled network architecture ( Fig 6F). This suggests that the probability of adhesion stability and ultimately maturation is higher with bundled architectures relative to both crisscrossed and random distributions of actin filaments (Fig 6F).
Collectively, our results indicate that the architecture of the actin cytoskeleton does not modulate the amount of ligand binding and integrin clustering. However, actin network architecture determines the physical distribution of ligand-bound integrins in adhesions, with bundled actin filaments increasing the probability of adhesion stability, consistent with previous experimental observations [50].

Discussion
Since our experiments show that different integrin types exist in nascent and mature adhesions (Fig 1B), a computational model is here developed in order to understand if differences in nanoscale physical properties of integrins reflect on adhesions. This is largely untested by experimental approaches because it is very challenging to simultaneously distinguish between integrin types and isolate their nanoscale physical properties. The model is used to study how ligand binding affinity, rate of integrin activation, strength of pairwise interactions, bond kinetics, as well as the architecture of a network of actin filaments modulate integrin organization in adhesions and stress transmission. Our results collectively show that ligand binding and integrin clustering are competing mechanisms and that bundled actin networks favor adhesions stability, and ultimately maturation.
The model is developed through three consecutive stages of increasing complexity: (i) simulations of single-point integrins diffusing on a quasi-2D surface and switching between active and inactive states, binding ligands, and interacting laterally; (ii) incorporation of an implicit actin flow and integrin/ligand catch bonds kinetics; (iii) binding of integrins to semi-flexible actin filaments in either random, bundled, or crisscrossed architectures. At all stages, we distinguish between β-1 and β-3 integrins, by using either exact, experimentally detected physical parameters, realistic fold differences between the two, or estimates from previous free energy calculations.
For high E IL , many active integrins bind ligands and the fraction of integrins that can diffuse, and cluster, is reduced (Fig 3D-3E). Accordingly, this happens when the fraction of β-1 integrins is higher than that of β-3 integrins (Fig 5A), since β-1 integrins have higher ligand- binding affinity than β-3 integrins. By contrast, with many free diffusing integrins that have low E IL , and are less likely to bind ligands, the fraction of integrins that can encounter each other, and cluster is enhanced and reduces ligand binding (Fig 3D-3E). This happens when the fraction of β-3 integrins is higher than that of β-1 integrins (Fig 5A) and also results from the higher diffusion coefficient of β-3 integrins with respect to β-1 integrins [43,44]. The result that ligand binding and integrin clustering are competing mechanisms is consistent with a kinetic Monte Carlo model showing that thermodynamics of ligand binding and dynamics of integrin clustering interplay [46]. Our model reproduces this competing process over the same range of ligand binding affinities and strength of pairwise interactions.
The molecular mechanisms resulting in integrin lateral clustering remain controversial. However, several lines of evidence have suggested that β -3 integrins assemble clusters more easily than β-1 integrins. For example, activation of β -3 integrins induces formation of clusters with recruitment of talin [7], while β-1 integrins require recruitment of many more signaling components in order to form clusters, such as FAK [57]. β-3 integrins cluster in response to talin binding without a concomitant increase in affinity [24], while β-1 integrins cluster only when extended [58]. Moreover, previous studies in U2OS cells showed that β-3 integrins cluster on both β-3 and β-1 integrin ligands, while β-1 integrin clusters are present in adhesions only on β-1 ligands [59]. Our result that β-1 integrins are correlated with ligand-binding while β-3 integrins are mostly responsible for clustering is consistent with the observation that clusters of β-1 integrins are present only on β-1 ligands, possibly because, in this case, ligand binding and not pairwise interactions facilitate adhesions assembly. Further evidence that β-3 integrins assemble more easily than β-1 integrins is provided by the reported spatial and functional segregation of the two integrin types. β-1 integrins translocate from the cell periphery to the cell center to withstand higher tensions, whereas β -3 integrins remain at the cell edges to do mechanosensitive activities via dynamic breakage and formation of multiple bonds with the substrate and with one another [38].
Our model also shows that the number of integrins per cluster, computed as the fraction of clustered versus ligand bound integrins, is in the range of 2-15 particles, depending on E II and E IL (Fig 3F). This value is comparable to the experimentally estimated number of integrins in nascent adhesions, between 5-7 [10]. By varying ligand density in the model, the ratio between clustered and ligand bound integrin particles does not vary significantly (S1 Fig), suggesting that the average number of integrins per cluster in nascent adhesions is not modulated by ligand concentration, consistent with previous experimental observations [10].
In the presence of actin flow, the fraction of β-1 integrins is positively correlated with ligand binding (Fig 5A). Above a threshold actin flow, however, ligand binding is almost suppressed, independent from relative amounts of β-1 and β -3 integrins (Fig 5A), because of faster ligand unbinding from both integrins. This reduction in bound ligands corresponds to a drop in the average tension per integrin upon increasing actin flow (Fig 5B). The biphasic response of tension to actin flow was previously observed experimentally [60] and is consistent with models of adhesion clutch assembly and rigidity sensing [61]. By increasing the fraction of β-1 integrins, a reduction of lateral integrin spacing is observed with our model (Fig 5C). Previous studies on the lateral separation of integrins in adhesions reported that a minimum spacing of 70 nm is required to form stable adhesions [62]. This value corresponds in the model to a minimum of 20% β-1 integrins (Fig 5B). This value represents a prediction from our CG model that can be experimentally tested in the future. When only β-3 integrins are used in the model, their lateral separation, upon binding ligands, is about 120 nm (Fig 5C), supporting the notion that β-1 integrins are needed to form stable adhesions. This is consistent with the fast binding/ unbinding dynamics of β-3 integrins previously observed in experiments [28].
By incorporating a positive feedback between actin filament engagement and integrin activation, as observed in [63,64], the competition between clustering and ligand binding is maintained in all actin architectures (Fig 6D-6E). This positive feedback represents the functional link between cytoskeleton and adhesions, where an increase in the probability of ligand binding results from binding actin, via an inside-out pathway [12,65]. Our model shows that the number of ligand-bound integrins with an average separation below 70 nm is enhanced with a bundled architecture (Fig 6F), suggesting that this configuration favors adhesion stability, and ultimately maturation [50]. When integrins bind a bundled network, they are likely to re-bind in close proximity because bundled filament architectures present filaments that are spatially closer than filaments of crisscrossed or random networks, forming a spatial trap for the receptors. Of interest for future studies is mimicking conditions of actin filament turnover, in order to understand how a dynamic cytoskeleton can interplay with integrin mixing in forming nascent adhesions. This will help understanding outside-in pathways, where, for example, adhesions formation modulates actin filaments polymerization. A further extension of the model will incorporate dynamic ligands, interconnected by a fibrous extracellular matrix that deform under tension. We will study how adhesions formation can change ligand localization and how this, in turns, affects adhesions morphology.
Previous computational studies of integrin dynamics range from all-atom simulations and enhanced sampling methods for understanding integrin activation at the level of individual molecules [66][67][68], to lower resolution coarse-grained [61,[69][70][71][72][73][74], lattice-based [75,76], diffusion-reaction algorithms [77] and theoretical models [78] for multiple integrins. With respect to the previous lower resolution models of multiple integrins, our new model allows us to directly incorporate properties of different integrin types, as detected experimentally (Fig 1B). The particle-based implementation scheme of our model is similar to that of other software for modeling the cytoskeleton, such as Cytosim [79] and Medyan [80]. However, important differences exist. In contrast to Cytosim, an explicit implementation scheme is used here because our time step, combined with the limited number of simulated particles (a few hundreds), allows us to achieve time scales of a few minutes, that are relevant for adhesion assembly, without excessive computational cost. In addition, in contrast to Medyan, our model does not have a scheme for solving stochastic reaction-diffusion equations, but instead focuses on the mechanics of particle interactions and displacements under deterministic and Brownian forces.
To conclude, with our highly coarse-grained model based on Brownian Dynamics, we extend the scope of previous theoretical and computational studies of integrin-based adhesions formation, by testing how differences in nanoscale properties of β-1 and β-3 integrins impact ligand binding, clustering and transmission of traction stress. By coupling physical parameters (such as diffusivity) together with chemical (i.e., affinity and receptor pairwise interactions) and mechanical (bond kinetics) parameters, and by using an explicit actin cytoskeleton, our model shows that nascent adhesions assembly can be finely tuned by differences in nanoscale physical properties of integrins. The CG model ultimately demonstrates that nanoscale differences in integrin dynamics are sufficient to determine ligand binding and integrin clustering. By incorporating dynamics of individual integrins in an explicit way, our model provides results that are consistent with a number of previous independent experimental observations, revealing important insight into the molecular origins of adhesion organization and mechanics. Taken together, our modeling results support the general view that a cell can control integrin expression to determine morphological and dynamic properties of adhesions.

Methods
In order to characterize how nanoscale physical properties of integrins impact the assembly of nascent adhesions, we developed a highly coarse-grained computational model based on Brownian Dynamics. The model is agent-based in the way sequential dependencies regulate interactions between integrin/ligand, integrin/integrin, and integrin/actin. Integrins, ligands, and actin filaments are explicit particles, while the cell membrane is implicit. Solvent molecules mimicking cytoplasmic effects are replaced by stochastic forces, depending on cytoplasmic viscosity. Inactive integrins diffuse and, when active, can bind ligands and interact laterally. When integrins are bound to ligands, they can engage actin filaments. The interaction between integrins and actin filaments locally increases integrin activation rate, ultimately resulting in a positive feedback between actin binding and ligand binding [61,81].
In order to distinguish between β-1 and β-3 integrins, we examine the effect of varying integrin activation rates, motility, ligand binding affinity, clustering, and bond kinetics (Fig  2A-2C). By varying the relative amounts of β-1 and β-3 integrins, we analyze fractions of ligand-bound integrins, clustered integrins, and average tension on integrin-ligand bonds (Figs 3-5). Moreover, we study the effect of different actin filaments architecture on adhesions morphology (Fig 6).
The model is an extension of our mechanosensing model [5] but differs from it in several ways. First, each integrin exists in either active or inactive state, determined by activation and deactivation rates. Second, the model incorporates tunable parameters for integrin physical properties, allowing us to discriminate between integrin types. Third, explicit semiflexible actin filaments are included.

Computational domain
The computational domain includes two systems: a square bottom surface, of 1 μm per side, and a rectangular 3D domain above the surface, with dimensions 1 x 1 x 0.04 μm (Fig 2A). The bottom surface mimics the substrate; the lower side of the rectangular domain mimics the ventral cell membrane above the substrate, while its inside space represents a 40 nm thick cytoplasmic region where actin filaments diffuse beyond the ventral membrane (Fig 2A). The cell membrane is separated from the substrate by 20 nm, a dimension characteristic of active integrin headpiece extension (Fig 1A) [82]. Within the cell membrane, integrins diffuse in quasi-2D and are restrained in the vertical direction by a weak harmonic potential with spring constant 100 pN/μm, mimicking membrane vertical friction. In the cytoplasmic region, a repulsive boundary is used on the top surface, to avoid filaments crossing the boundary. Periodic boundary conditions are applied on all lateral sides of the domain, in order to avoid finite size effects.

Integrin and ligands representation
The model considers a given number of ligands on the substrate, randomly distributed and fixed in space. We use a ligand density of 1000#/μm 2 , of the same order of that used in a previous model of adhesions assembly [72]. Integrin density on the cell membrane is~100#/μm 2 [5]. Integrins are single-point particles, that are initially randomly distributed and diffuse over the course of the simulations. Integrin diffusion coefficient is D = 0.1 μm 2 /s for β-1 integrins and D = 0.3 μm 2 /s for β-3 integrins [43]. Introducing volume exclusion effects between integrins, in the form of a weak repulsion between nearby particles (1 pN force), does not change the fraction of ligand-bound integrins, their average separation, the mean tension per integrin and its distribution (see S2A-

Actin filament representation
Semiflexible polymers represent actin filaments as spherical particles connected by harmonic interactions. Filaments have fixed length of 0.5 μm, corresponding to 6 beads separated by 0.1 μm equilibrium distance. The model of actin filaments is explained in detail in [83]. Actin filament beads are subjected to both stochastic and deterministic forces. Stochastic forces on the i-th bead are random in direction and magnitude in order to mimic thermal fluctuations and satisfy the fluctuation-dissipation theorem: withÎ a;b being the second-order unit tensor [84] and μ being a friction coefficient equal in three directions. Deterministic contributions come from bending and extensional forces on the filament beads. The bending force is computed as: where l p = 10 μm is actin filament the persistence length, N is the number of beads in a filament (N = 6) and t i ¼ ðr jþ1 À r j Þ jr jþ1 À r j j . The extensional force on filaments beads is computed as: where l 0 is the equilibrium length of 0.1 μm, k is the spring constant of 100 pN/μm. Each spherical particle of a filament represents a binding site for integrin and each binding site can interact with multiple integrins.

Agent-based algorithm
In order to mimic hierarchical formation of nascent adhesions [10], the algorithm incorporates sequential interactions between integrins, ligands and actin filaments. First, we simulate a system composed of only integrins and ligands, in order to explore the ways in which integrins cluster and bind ligands in an actin-independent way. Then, we add actin filaments and study the effect of actin network architecture on adhesions formation.
Integrins switch between inactive and active states, with rates of activation and deactivation k a = 0.5 s -1 and k d = 0.0001 s -1 , of the same orders of those previously estimated [44,45,72]. Activation probability corresponds to P a = k a dt, with time-step dt = 0.0001 s, as the time of the smallest simulated phenomena. This large time-step is allowed because the extensional stiffness of actin filaments, 100 pN/μm, is smaller than the real actin filaments stiffness, of about 400 pN/nm [85]. Upon activation, integrins can interact with free ligands, using a harmonic potential (with equilibrium separation 20 nm and spring constant 1 pN/μm), and cluster with other active integrins, depending on relative distances. Ligand binding occurs within a threshold distance of 20 nm, which reflects the extension of the open conformation of α IIb β 3 integrin away from the membrane [82]. Each integrin can bind only one ligand, and each ligand can bind only one integrin, mimicking binding sites specificity. Clustering occurs below a threshold of 30 nm, a value of the same order of the integrin-to-integrin lateral separation observed experimentally [82] and one order of magnitude lower than the minimum separation between individual adhesions [86].
The probability of integrin deactivation is P d = k d dt. Once inactive, integrin loses its connections with ligands and other integrins. Integrins unbind ligands with dissociation probabilities depending on their affinities: P ¼ le E IL dt, using a prefactor λ = 1 s -1 , for simplicity. They break later connections with probabilities inversely proportional to strength of pairwise interaction: P ¼ le E II dt. For β-1 integrins, we use high affinity,~9 k B T; for β-3 integrins we use lower affinity, 3-5 k B T.
Ligand-bound integrins can establish harmonic interactions with semiflexible actin filaments below 5 nm, approximating the size of the intracellular integrin tails [82]. Since the exact molecular composition of the layer between integrin and actin can contain up to~150 different proteins [87], a detailed modeling representation is not possible. Therefore, interactions between integrin and actin are approximated by harmonic potentials with equilibrium distance of 3 nm and spring constant of 1 pN/μm. These interactions simplify the~40 nm layer of adhesion molecules, including vinculin, talin and α-actinin, and is consistent with the level of details of the simulations, where harmonic interactions are used to connect particles within 20-100 nm.

Brownian Dynamics simulations via the Langevin Equation
Displacements of integrin and actin filament particles are governed by the Langevin equation of motion in the limit of high friction, thus neglecting inertia: where r i is a position vector of the ith element, z i is a drag coefficient equal in three directions, t is time, F i is a deterministic force, and F T i is a stochastic force satisfying the fluctuation-dissipation theorem [88]. F i is the sum of forces resulting from interactions of integrins with a ligand and/or other particles in the system, and actin flow in a direction parallel to the substrate. Positions of the various elements are updated at every time step using explicit Euler integration scheme: Integrin/ligand unbinding formalism Since contraction forces are not needed for the assembly of nascent adhesions [5,10,89], our computational model only incorporates forces mimicking actin retrograde flow. In order to simulate actin flow and characterize distribution of traction stress at various flow rates, a constant force is applied on ligand bound integrins, along y ( Fig 2B). Lifetime of the bond between integrin and ligand follows the catch-bond formalism (Fig 2C), using: for β-1 integrins an unloaded affinity of 2 s and a maximum lifetime of 15 s; for β-3 integrins an unloaded affinity of 0.5 s and a maximum lifetime of 3 s. The parameters for the catch bond kinetics are from previous experimental characterizations [38,47,56]. Curves of bond lifetime versus tension are shown in Fig 2C. For β-1 integrins, we implemented an unbinding rate as a function of the force acting on the bond, F: For β-3 integrins, we used unbinding rate: The functional forms of the catch bonds were taken from a model that assumes a strengthening and a weakening pathway for the bond lifetimes, using a double exponential with exponents of opposite signs [90,91]. This model was also used for previous simulations of integrinbased adhesions [5].

Positive feedback between filament binding and integrin activation
To mimic promotion of integrin clustering upon ligand binding and actin filament engagement [81], we introduce a positive feedback between binding of integrin to a filament and integrin activation rate. In the model, integrins can bind a filament only if already bound to a ligand. Upon binding to actin, integrin activation rate is increased by 2 to 4% relative to its initial value. This assumption is motivated by recent evidence from TIRF experiments on T-cells, where it was demonstrated that actin binding and correct ligand positioning are needed for integrin activation [81]. The positive feedback between actin binding and integrin activation rate also represents conditions of inside-out signaling, with increased affinity for ligand binding induced by the cytoplasm [65]. We use the model with the positive feedback (schematics in Fig 2D) to test the effect of different actin architectures on ligand binding and clustering (Fig 6).
For bundled and crisscrossed actin filament architectures, we impose spatial restraints on filaments pairs. Bundled architectures have harmonic connections between beads of filament pairs that keep the filaments in parallel; crisscrossed architectures impose 90 deg angle between the axis of filaments pairs, that keep them almost perpendicular.
Cells were imaged using a 1.2 NA 60X Plan Apo water immersion lens on an inverted Nikon Ti-Eclipse microscope using an Andor Dragonfly spinning disk confocal system and a Zyla 4.2 sCMOS camera. The microscope was controlled using Andor's Fusion software.