Matrix feedback enables diverse higher-order patterning of the extracellular matrix

The higher-order patterning of extra-cellular matrix in normal and pathological tissues has profound consequences on tissue function. Whilst studies have documented both how fibroblasts create and maintain individual matrix fibers and how cell migration is altered by the fibers they interact with, a model unifying these two aspects of tissue organization is lacking. Here we use computational modelling to understand the effect of this interconnectivity between fibroblasts and matrix at the mesoscale level. We created a unique adaptation to the Vicsek flocking model to include feedback from a second layer representing the matrix, and use experimentation to parameterize our model and validate model-driven hypotheses. Our two-layer model demonstrates that feedback between fibroblasts and matrix increases matrix diversity creating higher-order patterns. The model can quantitatively recapitulate matrix patterns of tissues in vivo. Cells follow matrix fibers irrespective of when the matrix fibers were deposited, resulting in feedback with the matrix acting as temporal ‘memory’ to collective behaviour, which creates diversity in topology. We also establish conditions under which matrix can be remodelled from one pattern to another. Our model elucidates how simple rules defining fibroblast-matrix interactions are sufficient to generate complex tissue patterns.

The organization of extracellular matrix in the body is diverse and often becomes dysregulated in cancer. How such patterns emerge is poorly understood, in particular how the interplay between the fibroblasts, which produce the matrix fibers, and the fibers themselves affects emergent organization. We created a unique adaptation to the Vicsek flocking model to include feedback from a second layer representing the matrix. Fibroblasts were able to deposit, degrade and rearrange fibers and the fibers could also contribute to the vector of cell migration. We find that this feedback is sufficient to generate many complex tissue patterns, as seen in vivo. The model is able to encapsulate key fibroblast and matrix properties without compromising on scale, providing an insight into the co-evolution of two crucial components of tissues at the scale of millimeters.

Introduction
The extracellular matrix (ECM) provides support and enables the functioning of tissues within multicellular organisms [1,2]. Depending on the tissue, the ECM can be organized in a myriad of different ways [3][4][5] from highly aligned linear bundles in connective tissue such as tendons, to largely anisotropic meshwork in basement membranes [6]. Curved bundles of ECM fibers are found in mammary tissue that can be loose, diffuse curls or tightly coiled bundles, for example in the dermis [7]. The patterning of ECM has large consequences on tissue function and often becomes disregulated in disease. How these different topologies emerge and subsequently interconvert as organisms age and develop pathologies, such as cancer, remains largely unknown. ECM is produced by a variety of different cells, with fibroblasts playing a particularly prominent role [8][9][10][11]. Numerous studies have documented the biosynthetic pathways of ECM production and some of the principles of fiber assembly [12,13]. However, the relationship between cell behavior and the emergent patterns of matrix organization at the mesoscale (millimeter scale) is unclear. The problem is further complicated because cell behavior is instructed by the ECM [14,15], with cells aligning themselves with matrix fibers and responding to different matrix stiffness with altered cell signaling and cell migration [1,16,17]. Mechanoreciprocity and mechanotransduction have been shown to be crucial phenomena in determining tissue and tumor evolution [18,19]. The feedback between cell behavior and matrix organization therefore presents a 'chicken and egg' problem when trying to understand the development of higher order patterns. Does matrix structure exert a dominant role over cell behavior that then reinforces the matrix topology? Or does cell behavior initiate matrix patterning? The lack of matrix at the zygote stage of development suggests that cell behavior is somehow the initial determinant of matrix patterning, but it does not preclude an important role for feedback from the extracellular matrix in specifying pattern formation. Fibroblasts are largely responsible for producing, degrading and rearranging matrix components and enabling their assembly into larger complex structures. This latter process is highly dependent on the integrin family of matrix receptors that span the plasma membrane and connect the cell's cytoskeleton to the extra-cellular matrix. Many perturbations have been identified that lead to altered matrix anisotropy, including within integrins, integrin-associated molecules, and the actin cytoskeleton, generating simple patterns of alignment and disorder. There have been numerous studies documenting the diversity of fibroblast behavior [20,21] and the ability of cells to cooperate, often across many orders of magnitude [22][23][24][25][26]. Yet it remains unclear how complex tissue patterns develop and why particular perturbations alter such patterns. Tackling this problem experimentally is challenging due to the complexity of systems with many cells and matrix components, and the difficulty of seeing in real time how matrix is being organized. To understand how specific rules defining cell-matrix interactions can generate emergent behaviors, we therefore turn to computational modeling.
Modeling of such behavior requires the inclusion of discrete fibers, rendering continuum models less appropriate. A Cellular Potts model has been used in a previous study to explore cell migration in ECM [27]; however, it is difficult to model explicit fibers to behave in a physiologically reasonable way and to implement rules for cell adhesion specifically between the cells and individual fibers without making the system far more complicated. Simple models of nematically aligning self-propelled rods [28][29][30] serve as a good basis for modeling collective cell behavior but move in a continuum fluid and have a simplistic averaging for coordination between rods. Hybrid models with discrete fibroblasts and a continuum matrix have been used to explore scarring in wound healing [31,32] and fibroblast alignment [33]. Work by Dallon et al. [34] explores interactions between the cells and matrix. They report that matrix feedback can reduce overall alignment, which we confirm here. The authors also show that a strip of aligned fibers that is sufficiently thick can cause the rest of the matrix to align similarly. Our work builds on such hybrid models by incorporating cell-cell interactions and cell migratory noise explicitly. Further, we model fibers discretely, allowing for an understanding of how an individual fiber can influence cell motility.
The Vicsek model [35] has been used for many decades to explore collective behavior of animals at high density, most notably shoals of fish and flocks of birds. In this approach, migration of animals or cells is simulated with individual velocities. Crucially, the vector of migration is determined by the combination of the level of intrinsic persistence [36] of a cell (influenced by the magnitude of a Gaussian noise term [37]), henceforth called individual migratory noise, and the propensity for a cell to align its directionality with that of nearby individuals. Within certain regions of parameter space this can generate rapid and effective alignment. Several extensions of the Vicsek model and models with similar neighbor alignment [29,38,39] or velocity alignment [40][41][42] have been developed to tackle specific biological questions [43], for example, how Dictyoselium discoideum can aggregate as a result of chemotaxis towards pulses of cyclic adenosine monophosphate and when this switch to collective behaviour occurs [44]. However, in its simple form, the Vicsek model is not well suited to modelling the interplay between cells and ECM. Matrix persists over long time scales and can be both synthesized and degraded, whereas the alignment term in the Vicsek model only considers the vector of individuals spatially local at the current time step. To circumvent this, we have developed a two-layer adaptation of the Vicsek model. One layer represents migratory cells and is similar to the conventional Vicsek model, the second represents the matrix and records the cumulative effect of cells on matrix synthesis and degradation. This second layer contributes an additional term in the behavior of the cells that represents matrix feedback on cell behavior. By exploiting this model in combination with experimentation we explore the role of matrix feedback in the generation of higher order patterns of ECM.

Diverse matrix patterns are found in vivo and can be quantified
Second harmonic imaging of collagen fibers in different mouse tissues in vivo revealed various patterns of matrix organization (Fig 1a). The position of fibroblasts was determined by using mice with transgenic expression of a fluorescent nuclear marker, H2B-GFP, in fibroblasts. Of note, the long axis of the nucleus indicates the orientation of the fibroblast cell body and a clear correlation between fibroblast and matrix alignment is visible in stomach ECM and has been previously documented [5].
To quantify the matrix diversity observed in vivo, we derived five metrics of matrix organization: long-range alignment of fibers (LRA), short-range alignment (SRA) (S1a Collectively, these metrics can describe a range of matrix properties. For example, a swirl-like matrix has low long-range alignment whilst maintaining some shortrange alignment (dermis, Fig 1a). Such a matrix would also have high curvature. The percentage of high density matrix encapsulates how homogeneous a matrix is. If cells are corralled, hence channeled to cover specific regions, this would result in more high-density matrix (liver, Fig 1a). Fractal dimension measures the self-similarity of the matrix and the extent to which the one-dimensional fibers can fill two-dimensional space [45]. A diffuse matrix will have a higher fractal dimension than matrix with large gaps and a non-aligned matrix, with fibers oriented in many different directions will have a higher fractal dimension than an aligned matrix, for example the spleen as compared to the stomach in Fig 1a. Examples of these metrics are given in S1 Fig. Application of these metrics to both murine and human pathological matrix confirms the diversity of matrix organization. For easy visualization of the different metrics, we generated star plots with five axes, each reflecting a different metric (Fig 1b).

Model description
To explore the mechanisms by which fibroblasts generate matrix patterns and how divergent matrix topologies initially arise, we initially adapted the Vicsek flocking model with cells moving subject to individual migratory noise and influence by their local neighbors [35]. A fibroblast's individual migratory noise reflects its persistent migration as a result of cell polarization [46]. At time zero, cells were placed at random and each cell assigned a random orientation (Fig 2) and a constant speed drawn from a Gaussian distribution. The change in the orientation of a cell is computed as a weighted function of three terms: individual migratory noise, cell-cell guidance and matrix guidance. Unlike the original Vicsek model, cells were able to coordinate their behavior nematically [28] and the extent of coordination could be explicitly controlled by the cell-cell guidance term. Fibroblasts are modelled as a bead representing the head of the cell, followed by two beads of twice the diameter representing the cell body, followed by another smaller bead for the tail. This diamond-like shape has aspect ratio 1:3 and reflects a typical cell morphology [47]. For a cell i at time t, θ i (t) denotes its orientation, s i denotes its speed and (x i (t), y i (t)) denotes its position.
Individual migratory noise. Cells move along their long axis at a constant speed and the individual migratory noise of a cell is modelled as a persistent random walk θ i (t) + η i (t), where η i (t) is Gaussian distributed with mean zero. The positional X and Y-components of this term  can be defined as Cell-cell guidance. If a cell i is in direct contact with N cells, N > 0, then the X and Ycomponents of the effect of cell guidance on i are defined as whereỹ j ðtÞ ¼ ( y j ðtÞ if j y j ðtÞ À y i ðtÞ j ðmod pÞ < p 2 ; y j ðtÞ þ p 2 otherwise: This adaptation constitutes the ability of fibroblasts to align in a nematic manner. For greater tractability, an assumption of the model is that w c remains fixed for any non-zero value of N, so that, provided a cell is in contact with at least one other cell, the degree of flocking will be fixed. This means cell-cell guidance is independent of N for N > 0. The angle of cell-cell flocking is an average of all N of a cell's contacts. However, the confluence levels both in silico and in vivo mean that N is kept small in practice. Matrix guidance. If the head of a cell i is above a fiber k with orientation ϕ k , then the X and Y-component of the effect of matrix guidance of k on i are representing the ability of fibroblasts to move along fibers in a nematic manner. The X and Y-components of cell i are then written as a weighted function of the three terms: where w p = 1 − w c − w m , 0 � w p , w c , w m � 1. If a cell is not in contact with any other cell, or above any fiber then w c or w m respectively are set to zero. The orientation of cell i at time t + 1 is then computed as which must then be adjusted for quadrant of the arctan function so that Finally cell position is updated so that where v e is a proxy for volume exclusion defined by The distance between cells which is considered as an overlap is set by the user but for this work was set as the inner 75% of the cell's area. Matrix updates. We modified the Vicsek model to include a second layer, consisting of a grid of matrix fibers. This meant the model comprised two layers: a top layer of fibroblasts, and an underlying layer of extra-cellular matrix. The ECM layer was arranged in a grid and each grid point was associated with matrix 'fibers' that could be oriented in eight possible bins (in increments of π/8 radians within the range [0, π] to reflect the nematic behavior of the system). When the head of a cell moved over a grid point it would update the fibers contained within the grid point. When establishing the model, we considered where fiber deposition would occur. In S2 Fig we show experimentally that fibronectin is produced throughout the cell and initially deposited in the front half of the cell, well ahead of the nucleus. We therefore chose to have cells depositing matix fibers below their heads. Fibers could be generated, deleted or assigned a different directionality, representing the deposition, proteolytic degradation and rearrangement of the matrix, respectively (Fig 2).
Fiber deposition. A cell with orientation θ deposits fibers according to the deposition rate in the grid point below its head in the bin which θ falls in [5]. If π � θ � 2π, fibers are deposited in the bin in which θ − π falls, since fibers are apolar.
Fiber degradation. All bins in the grid point below a cell's head will be depleted according to the degradation rate.
Fiber rearrangement. Fibers from the two neighboring bins will be moved to the bin which θ falls in at the specified rearrangement rate. The bins wrap around, so that if θ falls in bin 1, rearrangement will happen with its two neighboring bins: bin 2 and bin 8.
Choosing fibers for matrix guidance. Our main innovation to the original Vicsek model is in adding the ability of cells to be guided by the matrix fibers below it whilst simultaneously producing and reorganizing these fibers. We hypothesized that the migratory behavior of fibroblasts would be influenced by the emergent matrix organization and that this feedback would influence the matrix topology that ultimately emerges. To investigate directly if fibroblasts were indeed guided by the underlying ECM, we analyzed the migration of fibroblasts plated on matrices with differing degrees of alignment. Specifically, aligning and non-aligning fibroblasts generated respectively aligned and non-aligned matrix over the course of seven days. These original fibroblasts were then removed and new aligning fibroblasts were randomly seeded on both the aligned and non-aligned matrix (S7 Text). We found that the fibroblasts moving on the aligned matrix had significantly increased persistence as compared with the same fibroblast phenotype moving on the non-aligned matrix (Fig 2a, p = 0.0064, twotailed t-test). This implies that it is not just the presence of matrix but the topology of the matrix that determines cell motility.
Having established experimental justification for matrix guidance we added this mechanism into the model (Fig 2b and 2c). A cell will flock with fibers in the grid point below its head, the bin is determined according to a Gillespie algorithm that takes into account the density of fibers in each bin; therefore, bins with the highest number of fibers have the highest probability of being chosen. The cell will then flock with the left bound of that bin.
Sequence of model updates. The order of events in the simulations at each time-step is given as: Model implementation and parameter choice is described in the Supporting Text.

An overview of quantification of matrix metrics
Here we briefly describe how the five metrics are determined. A full description of their derivation is given in S2 Text.
Long-range alignment (LRA) and short-range alignment (SRA). The alignment for an individual fiber is computed as the median angle of deviation between that fiber and all fibers within a pre-defined neighbourhood length; a longer length for LRA and a shorter length for SRA. The alignment of a matrix image at this specified neighbourhood length is given as the mean of the alignment of individual fibers [5].
High density matrix (HDM). Matrix images are thresholded against a pre-defined value. Pixels above this threshold are considered to be fibers, whilst the remaining pixels are background. HDM is then computed as the percentage of fibers in the image.
Curvature (Curv). For each matrix image, by using a Gaussian filter with a large radius, a mask is derived showing the principal matrix fiber structure. The mask lines are then traversed by increments of a pre-defined curvature window. Curvature at that point, is the change in angle between one line segment and the next. The Curv metric is then the average of curvature values for the entire image [48].
Fractal dimension (Frac). For each matrix image, by using a Gaussian filter with a small radius, a mask is derived showing more in detail matrix fiber structure than the curvature masks. We then computed the box counting fractal dimension of the mask. For boxes of a given size, the number of these boxes required to cover the mask is computed. For boxes with side-length of ever-decreasing sizes, the number of boxes required to cover the mask is computed. Box dimension is then defined as the limit of how the number of required boxes to cover the mask scales with box side-length [45,48].

Individual migratory persistence and cell-cell guidance alone can generate alignment
We initially explored the interplay between cells' individual migratory noise and cell-cell guidance (Fig 3a, S3 Fig). The fibroblasts still deposited matrix into the second layer of the model, but w m was set to zero, reflecting the cells moving completely independently of the fibers. The model enables us to see different fibroblast and matrix patterning emerging over time (S1-S4 Videos). The top row of Fig 3a shows that in the absence of cell-cell guidance only isotropic matrix is generated, regardless of the level of individual migratory noise. Increasing the cellcell guidance term leads to the formation of both anisotropic matrix and a spatially non-uniform distribution of both matrix and cells (Fig 3a, S3 Fig).
Crucially this analysis showed that individual migratory noise (η) and cell-cell guidance (w c ) alone were not able to generate the wide diversity of patterns seen in vivo (Fig 1) [6,7]. Interestingly, there was surprisingly little variation in curvature amongst patterns and longrange alignment (LRA) and short-range alignment (SRA) are highly correlated (Fig 3b and 3c,  S3c Fig). This analysis indicates that noise and cell-cell guidance can generate a limited diversity of patterns, with the main emergent pattern besides isotropic matrix being uniform alignment.
In order to test experimentally if matrix feedback was dispensable for alignment, as the model suggested, we treated fibroblasts that generate an aligned matrix with the 'functional upstream domain' (FUD) of Streptococcus pyogenes. FUD prevents effective matrix assembly from soluble fibronectin into insoluble fibrils [49]. Cell migration is instructed by this fibril form and not the soluble form. In this way treatment with FUD precludes the cells from being guided by the matrix and we are able to eliminate matrix feedback (Fig 3d). The addition of FUD to fibroblast cultures efficiently prevented matrix fiber bundling but did not alter the migratory persistence of fibroblasts (Fig 3e) or prevent the progressive alignment of fibroblasts over time (Fig 3f). Taken together, these analyses indicate that a limited repertoire of matrix patterns can be generated in the absence of matrix feedback.

Matrix feedback generates diverse matrix patterns
We next considered the effect of introducing varying levels of matrix feedback whilst keeping the level of cell-cell guidance fixed at a level that is consistent with previous work and reflects experimental results [5] (w c = 0.03 as indicated by the pink strip in Fig 3a). This value of cellcell guidance provides a physiologically plausible level of coordination between cells such that cells with no individual migratory noise will align and cells with high individual migratory noise will not align. Given the restricted selection of patterns that can be generated by cell-cell guidance alone, we fixed cell-cell guidance at a biologically justifiable level and varied matrix feedback [5]. Starting with zero matrix feedback, cells with low individual migratory noise produced an aligned matrix characterized by high long-range alignment, high short-range alignment and low curvature (Fig 4a and 4b, red box). Increasing matrix feedback for these cells generated a matrix with high short-range alignment but low long-range alignment, medium curvature and a high proportion of high-density matrix (Fig 4a and 4b, yellow box). As documented in S5 Fig, cells with high individual migratory noise produced a diffuse isotropic matrix with high Frac and Curv and low SRA and LRA, and a low HDM (Fig 4a and 4b, blue box). Increasing matrix feedback for these cells produced a swirl-like but diffuse matrix with high short-range alignment, low long-range alignment and mid-range percentage of high-density matrix and curvature (Fig 4a and 4b, green box). To investigate further the effect of matrix feedback we performed a comprehensive analysis varying individual migratory noise, cell-cell guidance and matrix feedback together. To obtain an overview of this exploration of parameter space, we used principal component analysis (PCA) to reduce the dimensionality of our matrix metrics [50] (Fig 4c). Dimensionality reduction into two principal components could explain 85% of the variance. The loadings plot shows that the first principal component is related to an increase in short-range alignment, high-density matrix and a decrease in fractal dimension, largely describing the effects of increasing matrix feedback. The second principal component related to an increase in longrange alignment with a small reduction in curvature and high-density matrix, resembling the effects of increasing cell-cell guidance. Importantly, the addition of the matrix feedback term expanded the diversity of matrix outputs visible on the PCA plot (note the expanded region covered by the black dots in Fig 4c). Pairwise analysis of the metrics revealed that the addition of matrix feedback could produce patterns in new areas of metric-space (S5a and S5b Fig). In particular, matrix could now be produced with low long-range alignment but high short-range alignment, high percentage of high-density matrix (30%) and a wide range of curvatures. This describes the spectrum of 'swirl-like' patterns created by matrix feedback.
Interestingly, matrix feedback could act as a secondary mechanism for generating shortrange alignment of cells with high individual migratory noise (S5b Fig, blue line). However, for cells with low individual migratory noise, matrix feedback actually antagonizes the matrix alignment that would result from the cell-cell guidance term alone. S5b Fig also shows that starting with two cell phenotypes (with high (blue line) and low (yellow line) individual migratory persistence respectively), as matrix feedback increases, the patterning of the two cell-types converges.
Further, an increase in matrix feedback causes the curvature of cells with low migratory noise to increase, fractal dimension of both cell types to decrease and high-density matrix of Matrix feedback enables diverse higher-order patterning of the extracellular matrix both cell types to increase (S5b Fig). These changes reflect the increased 'corralling' of the cells as they follow and reinforce tracks of matrix fibers that were laid down at early time points (S5 Video). Visually, the result is pronounced swirl-like patterning. The matrix serves as a memory component to the flocking model, causing cells that would eventually normally align through cell-cell guidance to feedback with non-aligned matrix that had been deposited early on, thereby reducing overall alignment.
Similarly, for simulations run at sub-confluence, matrix feedback enhanced diversity of patterns (S5c and S5d Fig). While it is well-documented that cell shape is dynamic through space and time [51,52], fibroblasts often display an elongated morphology [5,47]. We hypothesized that different cell shapes in our coarse-grained model would not cause significant perturbations to the system. We confirmed this by running simulations with three distinct cell shapes (S6 Fig) [47]. We excluded the possibility that curvature was caused by the number of grid points comprising the matrix or the number of bins at each grid point by running a subset of simulations with finer grid points and more bins (S2 Text, S7 Fig). In conclusion, these analyses show that matrix feedback increases the diversity of emergent matrix organization, and enables more accurate recapitulation of matrix observed in real tissue.

How fibroblasts organize the matrix affects pattern formation
We next explored the effect of varying the different components of the matrix organization. In the previous simulations, fibroblasts had just deposited a single fiber at each time step, which had simply increased the number of underlying matrix fibers oriented in the direction of cell migration. We now varied the effect of cell migration on matrix fibers in three different ways: cells could increase the number of fiber elements, cells could re-align existing fiber elements, and cells could degrade existing filament elements, mimicking proteolysis. This manipulation of the second layer of the model allowed us to observe how dynamic remodeling of the ECM might affect matrix patterns by mapping new points in metric space onto the original PCA given in Fig 4c. To reflect the idea that matrix feedback is likely to have a lesser role at early time-points when little ECM has been deposited, we modified the model to take into account matrix feedback as a function of fiber density (see Methods). The original points from An increase in deposition rate results in higher HDM and higher curvature, under both high and low feedback regimes. Further, higher deposition rate causes lower LRA as the matrix becomes so dense that there is a loss of order. A higher rearrangement rate leads to fibers being organized into thicker bundles, leading to more high-density matrix especially for cells with high migratory noise (compare yellow and orange boxes S8b and S9b Figs). The effects of matrix reorganization remain largely constant regardless of level of matrix feedback (S8 and S9 Figs). These data demonstrate that varying matrix deposition and re-organization further enhance the diversity of matrix patterns. However, the overall effects of varying how cells reorganize the matrix are subtler than the absolute level of matrix feedback.

Diverse matrix patterns found in vivo can be mimicked in silico with matrix feedback
We next sought to relate insights from our model to physiologically observed ECM in Fig 1, shown again in Fig 5a and 5b. Through our exploration of parameter space we inferred parameters that would be able to generate similar matrix patterns, using a comparable number of fibroblasts (N = 200 cells, corresponding to 10% confluence) to the in vivo systems (Table 1, S5 Text). By reducing cell speed and making the fiber grid twice as fine to increase precision of matrix fiber positioning, the model was able to recapitulate all of these in vivo matrix patterns across the same mesoscale (Fig 5c). Nevertheless, given that the in vivo matrix patterns are more prone to noise within their biological component features as compared to the more exacting in silico simulation studies, great care was taken to pace comparisons on an equal footing, see S6 Text for an in-depth discussion. An interesting discussion on the effects of cell speed on emergent patterning due to distance travelled by cells between samplings can be found in [34].
The corralled and spatially heterogeneous patterns observed in the dermis and liver suggest that matrix feedback is required to generate these patterns and their high curvature is indicative of high migratory noise. We have shown that low individual migratory noise and nonzero cell-cell guidance alone are sufficient to generate a highly aligned matrix (Fig 2) like the stomach. The individual strand-like pattern of the spleen suggests a very low level of cell-cell coordination. The starplots characterizing the in silico matrix match the in vivo starplots closely after normalization (Fig 5d). Importantly, the 'swirl-like' patterns seen in the dermis and liver required matrix feedback (Table 1).

Interconversion between matrix patterns is inhibited by high matrix feedback
Matrix patterns are known to change in cancer and in ageing [7]. Having studied the emergent generation of matrix patterns, we next sought to study the interconversion from one matrix pattern to another. Matrix organization revealed by Gomori trichrome staining of collagen in a normal region of human dermis is characterized by curved matrix fibers, similar to the mouse (Fig 6a). In contrast, a region of dermis from the same individual that is adjacent to a melanoma exhibits altered matrix organization (Fig 6a), with fibers that are more aligned near the melanoma (Fig 6b). Consistent with a previous report that targeted BRAF inhibition in melanoma activates stromal fibroblasts [53], a comparison of ECM organization in the same patient prior to therapy and post-therapy indicated that when the treatment is failing, matrix alignment increases and curvature and fractal dimension decrease (S10 Fig). These analyses show that matrix can be dynamically remodeled over time [6].
Having observed matrix interconversion in pathological samples, we considered it in our model. At high levels, matrix feedback would simply constrain fibroblasts to migrate along existing ECM and preclude changes in matrix pattern. Fig 2 demonstrated that fibroblasts exhibit matrix feedback, but we wanted to determine if this was permissive for matrix remodeling or might be so high that it prevented matrix remodeling. We turned to the model to establish the conditions under which matrix conversion could be achieved. To mimic the situation in melanoma, we explored if an in silico generated 'dermal matrix' (Fig 5) would be altered by changing the properties of the fibroblasts once the initial dermal matrix had been deposited. As the matrix at the melanoma border exhibited higher alignment (LRA and SRA) than the normal dermis, we reduced the noise of the fibroblasts to zero in the model following the phase of initial matrix deposition. In line with expectation, the model predicts that high levels of matrix feedback will lead to fibroblast migration being channeled into following existing matrix, thereby preventing matrix remodelling (Fig 6c, blue line). However, lowering matrix feedback to 0 permitted a gradual transition from 'dermal' matrix to a more aligned matrix (grey line). The emergent matrix showed hybrid features with retention of the original 'dermal' matrix (Fig 6c). We therefore explored how adding fiber rearrangement and degradation might interplay with matrix feedback and conversion between matrix patterns (Fig 6c,  dotted lines). At a middle level of matrix feedback (w m = 0.2, Fig 5c orange line), enabling fibroblasts to rearrange and degrade fibers produces a significantly more aligned matrix (p = 9e-04, two tailed t-test). These analyses demonstrate how, when transitioning between two different matrix types, cells can take instruction from the original matrix [34] and produce a hybrid matrix that can have different patterning to de novo matrix generation (Fig 6d). In particular, high matrix feedback precludes interconversion by limiting fibroblasts to following existing matrix structures. We therefore conclude that matrix feedback needs to be both nonzero but also not too high. It has the benefit of diversifying matrix pattern, but cannot be so strong as to prevent matrix remodeling.
To address this experimentally, we seeded aligning fibroblasts on a thick non-aligned matrix and non-aligning fibroblasts on an aligned matrix, and observed over four days how the cells behaved on the matrix, and the patterning of the new matrix they produced (Fig 7). The original matrix is shown in yellow, with the emergent organization of the fibroblasts at different times shown in phase-contrast below. The new matrix produced by these fibroblasts is shown in blue. In both cases the fibroblasts display an initial tendency to follow the original matrix (day 1), but then revert to their preferred phenotype (day 3-4), being able to ignore the original matrix and produce a new matrix on top. Together, our adapted two-layer nematic Vicsek model and supporting experimental analyses demonstrate that matrix feedback enables diverse emergent patterns, including curved matrix structure, and the strength of matrix feedback is not sufficient to lock the system indefinitely, therefore enabling transitions to occur over a timescale of days.

Discussion
Simple rules give rise to complex patterns in many natural systems, from Turing's work on how a zebra gets its stripes to fractal self-similarity in ferns. Previous work [5,23] has shown how cell-cell interactions can lead to alignment in vitro. However, fibroblast and ECM organization will co-evolve in vivo in both normal and pathological tissues. Whilst there have been in-depth studies on how fibroblasts organize matrix via contractile forces generating local alignment around a cell [54], or between two contractile cells [5], and the effects of matrix on the movement of individual fibroblasts, surprisingly little is known about how the interplay between these two fundamental mechanisms act in tandem to self-organize. The emergent organization of ECM has critical consequences for tissue function but a mechanistic understanding of how some of these more complex ECM patterns arise is lacking.
To address these issues, we developed a nematic two-layer Vicsek flocking model, explicitly modelling matrix fibers and cells. The extent of matrix feedback and rates of matrix fiber deposition, degradation and rearrangement could all be adjusted in the model. The model enabled us to see both layers, the fibroblasts and the matrix, evolving simultaneously, which remains highly challenging in experimental systems. The model predicts that varying just the cell-cell guidance component without any matrix feedback can generate aligned and isotropic matrix, a finding that we have confirmed experimentally [5]. Introducing matrix feedback increased diversity of tissue patterns, generating matrix with varying long-range alignment, short-range alignment, percentage of high-density matrix, curvature and fractal dimension in metricspace. Matrix guidance of cells occurs irrespective of when the matrix fibers were deposited, resulting in feedback with the matrix acting as temporal 'memory' to collective behavior which creates diversity in topology. This temporal dimension of the model together with the ability to read and write between the two layers of fibroblast and matrix are novel extensions to the Vicsek model. Importantly, many of the patterns generated by the model are highly reminiscent of in vivo tissues at the mesoscale across millimeters. Recent work on dense suspensions of active filaments has shown it is also possible to produce swirl-like patterns of active filaments for specific values of persistence length and Péclet number [55], however this is dependent on a high packing fraction. Conversely, our model is able to produce diverse patterning even at low confluence, closer to cell density in certain tissues in vivo.
It has been hypothesized that matrix feedback could cause fibroblasts to align [22]. Our experiments suggest that it is cell-cell interactions which drive initial alignment without a requirement for ECM, but that matrix feedback plays a crucial role in determining cell motility when the matrix is thick and can facilitate alignment under specific conditions. How matrix feedback affects cell speed is not explored here, however a change in cell speed could affect the cells interaction with matrix, affecting the memory component of matrix feedback. These more nuanced effects of matrix feedback would be a topic of future interest.
Whilst the relatively few rules of the system allows for mathematical tractability and modelling on the scale of millimeters, the model lacks some more subtle matrix features such as matrix cross-linking, which we expect would act to limit matrix reorganization, and environmental cues such as chemoattractants in wound healing [32,33] and in collective behaviour of dictyoselium [44]. Additional work in this direction could add further insights into the effects of these phenomena. It would be of particular interest to use this model to better understand the transition between matrix topologies as seen in development [10,11] and ageing. A natural extension of the model would be to move into three dimensions, in order to more clearly understand the in vivo context. In three dimensions cell motility is different [56,57], and we conjecture that matrix feedback would be even stronger due to cells moving through, rather than on top of, matrix. However, continuing to model at the mesoscale level in three-dimensions would be computationally expensive. Finally, it would be interesting to use a deep learning network to see how well the model matches to matrix evolution in vitro and in vivo and how the model parameters change in these environments. Unlike other modeling approaches which take the ECM to be a surrounding fluid [22, 28-30, 32, 34], our two-layer Vicsek model has met the need for explicit representation of cells and matrix fibers, whilst remaining computationally inexpensive.
By coupling cell-cell and cell-matrix interactions together, we have shown how simple rules can produce diverse tissue patterns. The model offers a potential hypothesis for how many of the complex patterns seen in vivo could arise at the mesoscale level. Our multi-layered flocking model unifies cell and matrix behaviors and provides new mechanistic understanding of the consequences of matrix feedback, including its importance for curved matrix patterns and how high levels of feedback hamper matrix remodeling. Such observations are likely to prove important as robust strategies in tissue engineering [2], where the role of ECM-mediated growth and repair is of central importance.

Software and visualization
The model is written in C++. Data analysis was carried out in R and image processing was done in Fiji (ImageJ). Visualization of the matrix is done using openGL. For each grid point, a line showing the orientation of the density of the most recently selected bin is shown, with a color gradation from white to black as a function of that density up to a cut off of 25, which corresponds to a black line. Example visual output is shown in S1 Fig.

Computing persistence
To avoid bias in computing cell persistence due to cell collisions, cells were tracked at low levels of confluence. The program Metamorph was used for cell tracking. A range of time intervals T was chosen over which persistence should be computed. At each point time t, the persistence was taken to be the mean directionality ratio of all cells, computed as

Varying matrix feedback as a function of fiber density
We computed the persistence of sub-confluent cells moving on glass and on a thick aligned matrix (Fig 2a) for windows of one and two hours. Cells tracked for at least one hour were recorded and a spline tracing the cell's trajectory was produced using the loess package in R (α = 0.5). This was in order to smooth the intracellular movement, which resulted in many small fluctuations in the trajectories. For each cell, the median persistence for each cell was computed as described above. Simulations of single cells with varying values of Var(η) were run and the persistence computed for one and two-hour windows. The experimental results were then matched with simulations using a least squares approach to select the most likely value of Var(η) given a cell's persistence. The distribution of noise in the simulations that matched the cells moving on glass was Var(η) = 0.13. Simulations were then run for the same level of noise and incremental values of w m between 0 and 1. The persistence of these simulations was then compared with the persistence of the cells moving on thick aligned matrix. Using a least squares approach, the most likely value of matrix feedback on a thick matrix is w m = 0.27. Fig 3e and 3f suggest that for the duration of the four-day FUD experiment (where there is zero matrix guidance) as compared to the control experiment, there was no difference in the persistence or the emergence of alignment. We therefore assume that the control fibroblasts undergo little to no matrix guidance with matrix fibers produced for these first four days. Defining this mathematically we say that for days 1-4, w m f(d) = 0, where d is fiber density.