A Theoretical Approach to Coupling the Epithelial-Mesenchymal Transition (EMT) to Extracellular Matrix (ECM) Stiffness via LOXL2

Simple Summary Epithelial-mesenchymal transition (EMT) is a key process in cancer progression through which cells weaken their cell-cell adhesion and gain mobility and invasive traits. Besides chemical signaling, recent studies have established the connection of EMT to mechanical microenvironment, such as the stiffness of extracellular matrix (ECM). LOXL2 is representative of a family of enzymes that promotes fiber cross-linking in ECM. With increased cross-linking comes increased stiffness, which induces EMT that can, in turn, elevate LOXL2 levels. As such, a positive feedback loop among EMT, LOXL2, and ECM stiffness can be formed. We built a mathematical model on a core biochemical reaction network featuring this feedback loop, and showed how strongly it drives EMT. We also illustrated mechanistically how cross-linking connects with stiffness, using a mechanical model of collagen (a major component of ECM). Using this theoretical framework, we demonstrated the heterogeneity of LOXL2/stiffness and its implications on migrating cancer cells that could seed metastasis, the growth of secondary malignant tumors. This framework can inspire experimental studies of more fine-grained mechanotransduction and biomechanical heterogeneity in cancers. Abstract The epithelial-mesenchymal transition (EMT) plays a critical role in cancer progression, being responsible in many cases for the onset of the metastatic cascade and being integral in the ability of cells to resist drug treatment. Most studies of EMT focus on its induction via chemical signals such as TGF-β or Notch ligands, but it has become increasingly clear that biomechanical features of the microenvironment such as extracellular matrix (ECM) stiffness can be equally important. Here, we introduce a coupled feedback loop connecting stiffness to the EMT transcription factor ZEB1, which acts via increasing the secretion of LOXL2 that leads to increased cross-linking of collagen fibers in the ECM. This increased cross-linking can effectively increase ECM stiffness and increase ZEB1 levels, thus setting a positive feedback loop between ZEB1 and ECM stiffness. To investigate the impact of this non-cell-autonomous effect, we introduce a computational approach capable of connecting LOXL2 concentration to increased stiffness and thereby to higher ZEB1 levels. Our results indicate that this positive feedback loop, once activated, can effectively lock the cells in a mesenchymal state. The spatial-temporal heterogeneity of the LOXL2 concentration and thus the mechanical stiffness also has direct implications for migrating cells that attempt to escape the primary tumor.


Introduction
Metastasis remains the most lethal aspect of cancer progression. Various steps of the metastatic cascade have been intensely investigated from a biochemical signaling perspective. This has led to the concept of phenotypic plasticity-the ability of cells to reversibly

A mechanical Model of Collagen Network
We model extra-cellular matrices that are mainly composed of collagen as a 2D diluted, partially phantom, triangular lattice network ( Figure 2A). The potential energy consists of terms that resists bond stretching and bending. Taking the spacing between nearest neighbor vertices as unity, then for the bond between -th and -th vertices: Two consecutive initially colinear bonds joined by a single vertex form a fiber segment, which has an elastic bending energy: where = 0 if bonds − , − are colinear. The total potential energy is = + . Typically, the bond bending coefficient is taken to be much smaller than the bond stretching one, due to the small cross-sectional area of the individual fibers.
Each vertex acts as a free hinge of cross-linked fibers. A non-phantom and non-diluted triangular lattice has a connection number of 6. That is, three fibers join at each vertex. Partial "phantomness" means that with probability , a random one of the three

A mechanical Model of Collagen Network
We model extra-cellular matrices that are mainly composed of collagen as a 2D diluted, partially phantom, triangular lattice network ( Figure 2A). The potential energy consists of terms that resists bond stretching and bending. Taking the spacing between nearest neighbor vertices as unity, then for the bond between i-th and j-th vertices: Two consecutive initially colinear bonds joined by a single vertex form a fiber segment, which has an elastic bending energy: where θ ijk = 0 if bonds i − j, j − k are colinear. The total potential energy is E total = E stretch + E bend . Typically, the bond bending coefficient is taken to be much smaller than the bond stretching one, due to the small cross-sectional area of the individual fibers. Each vertex acts as a free hinge of cross-linked fibers. A non-phantom and non-diluted triangular lattice has a connection number of 6. That is, three fibers join at each vertex. Partial "phantomness" means that with probability p phan , a random one of the three fibers detach from the vertex and no longer couple to the remaining two. The lattice network is also diluted, i.e., each bond is formed only with a probability p bond . As such, the total average connectivity for a vertex is: This network stiffens with increasing z, as demonstrated by the computed linear shear modulus depicted in Figure 2B. To calculate the linear shear modulus, we apply a small strain γ with Lee-Edwards periodic boundary condition (i.e., by horizontally shifting neighboring periodic images above and below the computational box), and thereafter relax the network by minimizing E total = E stretch + E bend . The linear shear modulus is then where A is the area per period. The same formula for G(γ) is also applicable in the nonlinear regime of large strains. fibers detach from the vertex and no longer couple to the remaining two. The lattice network is also diluted, i.e., each bond is formed only with a probability . As such, the total average connectivity for a vertex is: This network stiffens with increasing , as demonstrated by the computed linear shear modulus depicted in Figure 2B. To calculate the linear shear modulus, we apply a small strain with Lee-Edwards periodic boundary condition (i.e., by horizontally shifting neighboring periodic images above and below the computational box), and thereafter relax the network by minimizing = + . The linear shear modulus is then = , where is the area per period. The same formula for ( ) is also applicable in the nonlinear regime of large strains.

Modeling Local Crosslink and Stiffness Variation in Extra-Cellular Matrices
LOXL2 promotes crosslinks between fibers. We take this into account by assuming a direct correspondence between a local and the local LOXL2 concentration. This decrease in phantomization mimicking the increased cross-linking increases the connectivity and hence the stiffness. To measure the local stiffness in this lattice just described, we perform bead experiments in silico as in [41]. To measure the stiffness at some location, we embed a bead therein and apply external harmonic "optical traps" along four orthogonal directions, , . The extra potential energy due to the trap is: where ∥ and ∥ are the position of the bead and trap bottom respectively, measured with respect to the initial location of the bead along or axis whichever is parallel to the trap gradient. The total potential energy = + + is then minimized with respect to the bead and vertex degrees of freedom. The local compliance (inversely related to stiffness) is defined as: Such a measurement can be repeated at different locations on the same network, one bead at a time. Fixed boundary conditions are applied; using a periodic boundary as in the bulk shearing would result in the whole lattice translating with no energy penalty. To

Modeling Local Crosslink and Stiffness Variation in Extra-Cellular Matrices
LOXL2 promotes crosslinks between fibers. We take this into account by assuming a direct correspondence between a local p phan and the local LOXL2 concentration. This decrease in phantomization mimicking the increased cross-linking increases the connectivity z and hence the stiffness. To measure the local stiffness in this lattice just described, we perform bead experiments in silico as in [41]. To measure the stiffness at some location, we embed a bead therein and apply external harmonic "optical traps" along four orthogonal directions, ±x, ±y. The extra potential energy due to the trap is: where d bead and d trap are the position of the bead and trap bottom respectively, measured with respect to the initial location of the bead along x or y axis whichever is parallel to the trap gradient. The total potential energy E total = E stretch + E bend + E trep is then minimized with respect to the bead and vertex degrees of freedom. The local compliance (inversely related to stiffness) is defined as: Such a measurement can be repeated at different locations on the same network, one bead at a time. Fixed boundary conditions are applied; using a periodic boundary as in the bulk shearing would result in the whole lattice translating with no energy penalty. To avoid artifacts due to edge effects, the measurement is performed only over a slightly smaller inner region of the whole lattice.

Modeling Positive Mutual Feedback between Cells and ECM
We build on the previously established set of ordinary differential equations modeling the core EMT biochemical network in order to include the new feedback loop discussed in this text ( Figure 4A). We solve for the steady state under different levels of external signaling and draw the bifurcation diagram using the software MatCont (version matcont6p11, https://sourceforge.net/projects/matcont/ (accessed on 15 March 2020)) with a custom MATLAB script.

Simulating EMT/MET Dynamics of a Migrating Cell
We take the spatial distribution of LOXL2 concentration and local stiffness from other sections of this work. Small number of cells are simplified as migrating in a constant velocity. The set of ODEs describing core EMT biochemical kinetics still applies, but due to the cell motility, the mechanosensing signal is now cell position-dependent and thus timedependent. Then, the ODEs can be solved by common solvers such as MATLAB's ode45.

LOXL2 Correlates Positively with EMT Markers and a More Mesenchymal Signature
First, we investigated the correlation between the expression levels of LOXL2 and that of key EMT markers for tumor samples, using various TCGA datasets (BRCA, COAD, COADREAD, OV), and for cancer cell lines using the Cancer Cell Line Encyclopedia (CCLE) dataset. LOXL2 consistently correlated positively with ZEB1 but negatively with CDH1 (a standard epithelial marker), across all datasets ( Figure 1). We have performed additional analysis of correlating LOXL2 with CDH1 (epithelial marker) and ZEB1, SNAIL2, SNAIL1, LOX, VIM, CDH2 (all mesenchymal markers/inducers). Across the datasets shown in Figure 1, LOXL2 correlates negatively with CDH1 and positively with others (see Supplementary Materials Table S1).
Next, we calculated the EMT scores for these cell lines and tumor samples using three different EMT scoring methods (76GS, KS, MLR) and checked their correlation with LOXL2. The 76 GS method has no specific pre-defined range of values and a higher 76GS score associates with an epithelial phenotype. Unlike the 76 GS method, the MLR and KS methods have predefined scales for EMT scores. MLR and KS score EMT on a range of [0, 2] and [−1, 1] respectively, with higher scores indicating a mesenchymal phenotype. The 76GS score showed a negative correlation with LOXL2 and as expected, the other two scores (KS, MLR) showed a positive correlation with LOXL2 across all tumor datasets ( Figure S1). This analysis signifies that the LOXL2 gene expression levels are associated with a more mesenchymal phenotype in tumor datasets. However, these patterns were not apparent in the CCLE dataset.

Cross-Linking Effect of LOXL2 Stiffens the Collagen Network
Given this correlation, we were motivated to develop a modeling framework that could account for this effect given the known interactions between these molecular players. As already mentioned, it is well-known that LOXL2 can increase the stiffness of extracellular matrix by cross-linking collagen. In recent years, there has been considerable effort devoted to formulating computational models of fibrous matrices typical of what would be expected in the stromal regions surrounding tumors. One popular approach is to consider a diluted lattice representation of the mechanics (triangular for a 2D system, face-centered-cubic in 3D) with both stretching and bending energies governing deformations of each link [42,43]. An example of one such diluted triangular lattice in 2D, is shown in Figure 2A. As discussed at length in the literature, the mechanics of this model can be highly non-linear, exhibiting strong shear-stiffening [42]. Recently, Mackintosh and co-workers have extended this basic approach to take into account variations in cross-linking [44,45]. Specifically, they introduced the idea of phantom crossings, sites at which some fibers (i.e., the continuing links) that pass through that site do not interact with other fibers.
For all the aforementioned lattice mechanical models, there exists a special critical value of the average connectivity z (number of bonds connected) of each lattice point, the so-called isostatic point z IP , whose value equals two times the dimensionality. ECM that is under normal physiological conditions is considered to lie below the isostatic point. For example, in the 2D models by Mackintosh and co-workers, there is always full phantomization. i.e., every site is phantomized, such that the connectivity is always below 4 = 2 × 2. This then gives rise to the aforementioned strain-stiffening and explains similar findings from experiments on reconstituted collagen I networks [44,45].
We adopt the idea of phantomization but make it tunable, in order to capture the upregulated cross-linking in cancerous tissues. For simplicity, we do this using the 2D triangular lattice framework (see Methods) as the basic physics of the stress response is independent of spatial dimension [45]. If the probability of having a node with this assumed "phantom-like" behavior is p phan and the probability of a link being present is p bond , the average connectivity becomes: As already mentioned, when z is below the isostatic point z IP = 4, the network stiffens with sufficiently large strain by entering the nonlinear regime; this is the aforementioned physiological strain-stiffening. Above this point, the network exhibits large stiffness even in the small-strain linear regime (Figure 1 in [45] presents an illustrative diagram). Either by strain-stiffening or by increasing z above this point, the stiffness becomes dominated by the (large) stretching modulus and is roughly independent of the (small) bending energy term. In this work we focus on the increase in z by decreasing p phan (via LOXL2) enabling the linear regime to attain high stiffness, as what is likely in cancerous pathology. Note also that within the correct range of parameters, p bond can be varied to account for varying fiber density.
Using this framework, we model the molecular effect of LOXL2 as causing a decrease in the value of p phan . Figure 2B shows a set of response curves for the stiffness as we vary this parameter for different values of p bond . As can be seen, the stiffness exhibits a significant rise once we have exceeded the critical value z IP by sufficiently lowering p phan . For our purposes, this means that as LOXL2 concentration is increased past a threshold value, it will give rise to an increasingly stiff ECM. As the threshold clearly decreases with p bond , it will be easier to create stiff ECM at higher values of fiber density.

Effects of Spatial Distribution of LOXL2
To get some idea of what might happen in practice for a localized tumor, we imagine that cells in some spherical region of radius R form a tumor of a constant cell density ρ and are emitting LOXL2 at a constant rate α. These molecules then diffuse with diffusion constant D and also decay at rate λ. In steady state, this creates an exponentially decaying concentration once we leave the tumor mass itself, namely: where the decay rate is d = √ D/λ. By assuming that this concentration modulates p phan via: we can obtain a map of the factor by which the local compliance is decreased, corresponding to stiffness being increased, due to the cross-linking effect (See Section 2). One such map for nominal parameter values is depicted in Figure 3. It shows how ECM outside of the tumor boundary is stiffened, and demonstrates a potential non-cell-autonomous biomechanical effect on cells neighboring the bulk tumor. It would be useful in the future to measure this behavior directly by monitoring both LOXL2 concentration and compliance maps for an in vitro tumor spheroid.
we can obtain a map of the factor by which the local compliance is decreased, corresponding to stiffness being increased, due to the cross-linking effect (See Section 2). One such map for nominal parameter values is depicted in Figure 3. It shows how ECM outside of the tumor boundary is stiffened, and demonstrates a potential non-cell-autonomous biomechanical effect on cells neighboring the bulk tumor. It would be useful in the future to measure this behavior directly by monitoring both LOXL2 concentration and compliance maps for an in vitro tumor spheroid. Figure 3. Local compliance measurement predicted by the triangular lattice network model. Cells inside the inner circles secrete LOXL2, which is then subject to diffusion and degradation, and modifies local collagen network connectivity. The outer grey parallelogram signifies the measured region, which is some distance away from the true computational boundary to avoid edge effects. The network is not rectangular because it follows the base vectors of the triangular unit cell.

LOXL2-Mediated Feedback Loop Drives EMT
The ECM stiffness can couple directly to EMT phenotype via a variety of signaling pathways [25,28,46]. Elevated LOXL2 level increases ECM crosslinking and the increased mechanical stiffness is then input to SNAIL, as indicated by the mechano-sensing studies mentioned above [25,26]. Thus, we can expect that as the matrix stiffens, cells will assume more of a mesenchymal phenotype, thereby increasing the production levels of the ZEB1 transcription factor. Further, [30] shows that increasing ZEB1 will upregulate LOXL2. There is then a positive feedback loop which is expected to stabilize mesenchymal phenotypes. Note that this non-cell-autonomous loop is distinct from possible intracellular effects of increased LOXL2 levels.
To demonstrate this effect, we implemented this feedback loop in our previously validated model of the core circuit governing EMT that consists of two interconnected, double-negative feedback loops and exhibits epithelial, hybrid, mesenchymal states [47].
We formulate the regulation from ZEB to LOXL2 using Michaelis-Menten kinetics: Figure 3. Local compliance measurement predicted by the triangular lattice network model. Cells inside the inner circles secrete LOXL2, which is then subject to diffusion and degradation, and modifies local collagen network connectivity. The outer grey parallelogram signifies the measured region, which is some distance away from the true computational boundary to avoid edge effects. The network is not rectangular because it follows the base vectors of the triangular unit cell.

LOXL2-Mediated Feedback Loop Drives EMT
The ECM stiffness can couple directly to EMT phenotype via a variety of signaling pathways [25,28,46]. Elevated LOXL2 level increases ECM crosslinking and the increased mechanical stiffness is then input to SNAIL, as indicated by the mechano-sensing studies mentioned above [25,26]. Thus, we can expect that as the matrix stiffens, cells will assume more of a mesenchymal phenotype, thereby increasing the production levels of the ZEB1 transcription factor. Further, [30] shows that increasing ZEB1 will upregulate LOXL2. There is then a positive feedback loop which is expected to stabilize mesenchymal phenotypes. Note that this non-cell-autonomous loop is distinct from possible intracellular effects of increased LOXL2 levels.
To demonstrate this effect, we implemented this feedback loop in our previously validated model of the core circuit governing EMT that consists of two interconnected, double-negative feedback loops and exhibits epithelial, hybrid, mesenchymal states [47].
We formulate the regulation from ZEB to LOXL2 using Michaelis-Menten kinetics: . We neglect any delays in the ECM modification by LOXL2, and instead make p phan directly depend on LOXL2 levels. Since the stiffness is directly dependent upon p phan , the effect is equivalent to inputting a mechanosensing signal to SNAIL: we take this to be of the form: Note that although we can calculate the curve of ECM stiffness vs. p phan (Figure 2B), the precise quantitative form of mechanosensing is unknown and hence Equation (5) at present is not meant to be exact. This choice of a Hill function captures the fact that, for some intermediate p bond , p phan must decrease beyond the isostatic point before the mechanical lattice network stiffens. Qualitatively, Equation (5) combines Equation (3), Figure 2B, and the unknown quantitative form of mechano-sensing. The reaction network is schematically displayed in Figure 4A. Readers looking for the detailed parametrization are directed to the Supplementary Materials.
The threshold constant K LE in the Hill function (Equation (5)) determines the strength of mechanosensing feedback. Note that the absolute concentration of LOXL2 is irrelevant; only the ratio K LE c LOXL matters. As shown in Figure 2B, ECM of higher density stiffens more easily; it starts to exhibit high stiffness with a lower concentration of LOXL2. This is equivalent to a smaller threshold K LE . Interpreted differently, K LE can be viewed as the "measuring stick" for c LOXL . A smaller K LE amounts to a larger g LOXL k LOXL /K LE using this "measuring stick". Since the ordinary differential equations are valid both for a single cell or for a group of nearby cells of similar expression state, increasing LOXL2 could be attributed to higher secretion rates or higher cell density of the secreting cells. Specifically, if we assume that all nearby cells are in the same state, then the overall LOXL2 level scales with the cell density multiplied by the single cell secretion.
Shown in Figure 4B is a family of steady-state bifurcation curves of ZEB1 vs. external signal I, going from right to left with decreasing K LE . A few features are worth noting. First, the epithelial states with very low ZEB values are not altered in any significant manner; except for the very extreme case where the mechano-sensing signal is always at maximum regardless of LOXL2 concentration, and hence only the mesenchymal state survives (K LE = 0, see Supplementary Materials Figure S3). Next, the hybrid E/M states become increasingly rare as the feedback is ramped up. This is because having an intermediate value of ZEB is untenable if it can take part in the new positive feedback loop that we have identified. Finally, there is a large hysteretic effect whereby mesenchymal states can remain stable even as the inducing signal is turned off, as M cells bootstrap their own existence by increasing the ECM stiffness in their immediate vicinity. This type of dynamical effect has previously been suggested as a consequence of TGF-β secretion by mesenchymal cells [48,49] as well as a possible "epigenetic locking" of cells in a mesenchymal state [50,51], but the results here are more dramatic because of the very strong nonlinearity of the stress response.

Spatial-Temporal Correlation and Cancer Cell Migration
As discussed in previous sections, the secretion, diffusion and degradation of LOXL2 creates a spatially varying ECM; and the LOXL2-ECM positive feedback loop drives EMT. Upon attaining mesenchymal status, small numbers of cells may start to migrate away from the site of primary tumor and initiate the process of metastasis. Their speed is usually much larger than the rate of bulk tumor growth and ECM remodeling due to this growth [52][53][54]. Thus, one can picture those cells as migrating through a temporally steady but spatially heterogeneous medium that is able to supply a varying EMT-inducing signal to them and dynamically couple to their intracellular EMT circuits.

Spatial-Temporal Correlation and Cancer Cell Migration
As discussed in previous sections, the secretion, diffusion and degradation of LOXL2 creates a spatially varying ECM; and the LOXL2-ECM positive feedback loop drives EMT. Upon attaining mesenchymal status, small numbers of cells may start to migrate away from the site of primary tumor and initiate the process of metastasis. Their speed is usually much larger than the rate of bulk tumor growth and ECM remodeling due to this growth [52][53][54]. Thus, one can picture those cells as migrating through a temporally steady but spatially heterogeneous medium that is able to supply a varying EMT-inducing signal to them and dynamically couple to their intracellular EMT circuits.
More specifically, we imagine the ECM whose spherically symmetric LOXL2 distribution is described by Equation (2) and whose mechanical stiffness is mapped as in Figure  3. Cells are taken as traveling from the edge of the primary tumor, at a constant velocity in the radial direction, towards a nearby blood or lymph vessel. The relatively small number of these mobile cells presumably means that they do not significantly alter the ECM on their paths, nor is there significant growth of the primary tumor on this time scale, so the LOXL2-ECM feedback loop is negligible, contrary to the previous case of a cluster of cells in like states. In addition to some constant background biochemical signaling , the ECM supplies a mechanosensing signal as in Equation (5), but now varies as the cell migrates through the heterogeneous medium. The LOXL2 concentration is determined by the inherent traits of the primary tumor, and could be altered for example by LOXL2-suppressing drugs [55]. In Figure 5 we show the resulting EMT circuit dynamics inside such migrating cells for different values of . The cells maintain the mesenchymal states until a critical location at which sharp drop of ZEB indicates that the total external signal has fallen to the regime that supports only epithelial states. Arguably, any such cell must reach the blood or lymph circulation vessels before such a mesenchymalto-epithelial transition completes. Otherwise, it becomes stuck not far from the primary tumor, and most likely cannot seed the metastasis process-fortunate for the patient. More specifically, we imagine the ECM whose spherically symmetric LOXL2 distribution is described by Equation (2) and whose mechanical stiffness is mapped as in Figure 3. Cells are taken as traveling from the edge of the primary tumor, at a constant velocity in the radial direction, towards a nearby blood or lymph vessel. The relatively small number of these mobile cells presumably means that they do not significantly alter the ECM on their paths, nor is there significant growth of the primary tumor on this time scale, so the LOXL2-ECM feedback loop is negligible, contrary to the previous case of a cluster of cells in like states. In addition to some constant background biochemical signaling I, the ECM supplies a mechanosensing signal I ECM as in Equation (5), but now I ECM varies as the cell migrates through the heterogeneous medium. The LOXL2 concentration αρ is determined by the inherent traits of the primary tumor, and could be altered for example by LOXL2-suppressing drugs [55]. In Figure 5 we show the resulting EMT circuit dynamics inside such migrating cells for different values of αρ. The cells maintain the mesenchymal states until a critical location at which sharp drop of ZEB indicates that the total external signal has fallen to the regime that supports only epithelial states. Arguably, any such cell must reach the blood or lymph circulation vessels before such a mesenchymal-to-epithelial transition completes. Otherwise, it becomes stuck not far from the primary tumor, and most likely cannot seed the metastasis process-fortunate for the patient.

Discussion
Cancer metastasis is one of the most vicious traits of cancer. The ability to metastasize is related to the plasticity of cancer cells, of which the epithelial-mesenchymal transition (EMT) is a key axis. Recently, it has been shown that in several cancer types, there can

Discussion
Cancer metastasis is one of the most vicious traits of cancer. The ability to metastasize is related to the plasticity of cancer cells, of which the epithelial-mesenchymal transition (EMT) is a key axis. Recently, it has been shown that in several cancer types, there can exist a possible positive feedback loop along the SNAIL1-ZEB1-LOXL2 axis [30,[32][33][34][35][36][37][38][39]. Elevated levels of all three molecules are correlated with cancer malignancy [56][57][58][59][60]. Among them, LOXL2 is downstream to ZEB1 and has a known extra-cellular effect. It promotes crosslinking between collagen fibers, a major component of the extra-cellular matrix (ECM), thus stiffening and stabilizing the ECM. Increased stiffness can be fed back to SNAIL1 (i.e., increased stiffness can upregulate SNAIL levels) and further drives EMT.
In this article, we have examined the bioinformatics-based, biochemical, and mechanical aspects of this feedback loop. We do note that there can be additional intracellular effects of LOXL2 [37]. These are not taken into account in our circuit as they do not represent a significant change to the original core circuit with its assumed self-activation positive feedback loop of ZEB. We also do not take into account other mechanisms whereby LOXL2 (or other Lysyl Oxidases) can enhance tumor aggressiveness; these certainly include angiogenesis [61], pre-metastasis niche formation [62], drug resistance [63] and the growth of lymphatic vessels [64] at least partially due to upregulation of VEGF. These effects would need to be included in any attempt to comprehensively treat the role of LOXL2 in cancer progression and the overall efficacy of inhibitory drugs.
Studying gene expression patterns from TCGA and CCLE database revealed a positive correlation between LOXL2 and ZEB1, as well as a negative correlation between LOXL2 and E-cadherin. These findings are consistent with the notion that LOXL2 drives EMT and with the possibility that this driving takes place through the mechanical effect of increased collagen cross-linking.
To demonstrate this possibility, we adapted the triangular lattice network model typically used to study the mechanical properties of collagen networks (including local stiffness). Specifically, we associated bond dilution with collagen density, and vertex phantomization with collagen crosslinks. Under small strain in the linear regime, this network stiffens significantly only when the connectivity is increased past two times the dimensionality, in our 2D case, four. A simplified biophysical setup is simulated, where a geometrically localized cluster of cancer cells secretes LOXL2 that is then subject to diffusion and degradation. In the steady state, increased LOXL2 and stiffness can "seep" outside the range of the cluster. From the biological perspective, this is a non-cell-autonomous effect.
To have a more quantitative picture, we incorporated this interaction into the core EMT regulatory network previously developed [47]. The input of the mechanosensing signal is assumed to take the form of a high order Hill activation function. Different threshold constants in the Hill function thereby correspond to how much the secreted LOXL2 can stiffen the ECM. Under circumstances that can correspond to a very low collagen density, increased cross-linking by LOXL2 barely stiffens the network, the corresponding Hill function always stays below threshold, and the EMT bifurcation diagram is hardly affected. Under circumstances on the opposite extreme, the bifurcation is strongly shifted such that the purely mesenchymal states stabilize and dominate.
Another factor that can affect the previous result is the population structure of the cells. The fact that stiffness interacts with EMT provides a mechanism for cells to affect each other; a population of mesenchymal cells can push a nearby cell to also undergo EMT. Although it is in principle possible that a single cell could secrete enough LOX2 to alter its own microenvironment, it seems much more likely that this feedback sets in only in the presence of a sufficiently numerous population of cells. This hypothesis could easily be tested experimentally by varying cell density. This proposed mechanism implies that just as in the case of Notch signaling and of cell-cell paracrine coupling via release of TGF-β, the single-cell picture must give way to a more comprehensive tissue-scale analysis of EMT phenotypes and their spatial correlation [65,66].
Our modeling results offer an explanation of the experimentally observed correlations between the EMT-related transcriptional factors and extracellular enzymes. It also should inspire future studies integrating mechanics and cancer-related biology both experimentally and theoretically. Especially, we think that more accurate experiments mapping local matrix stiffness and ECM-modifying enzyme concentrations can test our model and can also open the door to a better understanding of the detailed roles of mechanical properties in cancer development. As an example, we also simulated the migration of individual mesenchymal cancer cells through such a heterogeneous, LOXL2-modified ECM. In this process, both spatial and temporal factors affect the "successful" seeding of metastasis, which depends on having cells navigate through a gradually diminishing EMT-inducing/maintaining stimulus, and attempt to remain mesenchymal in order to reach the blood or lymph circulation. More precise and detailed mapping of local stiffness, LOXL2 concentration, the effects of drugs, EMT status of cells as well as their time-dependence, will be key to deeper insights in these processes.
The current approach can be extended in a number of ways. For example, we focus on the cross-linking by LOXL2 while keeping the collagen deposition density as an external factor. However, there is evidence that cancer cells can also promote collagen deposition [30]. The diluted, partially phantom lattice network described in this work can encompass this additional effect via varying the bond dilution factor. We also could include the aforementioned intracellular role of LOXL2 by extending our baseline EMT circuit. One can also include more directly the effects of other members of the LOX family. These extensions will become useful once more quantitative data becomes available connecting the various pieces of this novel feedback mechanism. We hope that our predictions and demonstrated modeling capabilities will spur experimental groups to obtain key data regarding the mechanisms elucidated here.

Conflicts of Interest:
The authors declare no conflict of interest.