Theories of Mechanically Induced Tissue Differentiation and Adaptation in the Musculoskeletal System 1

It is well known that the mechanical environment affects biological tissues. The importance of theories and models that aim at explaining the role of the mechanical stimuli in process such as differentiation and adaptation of tissues is highlighted because if those theories can explain the tissue’s response to mechanical loading and to its environment, it becomes possible to predict the consequences of mechanical stimuli on growth, adaptation and aging of tissues. This review aims to present an overview of the various theories and models on tissue differentiation and adaptation, and their mathematical implementation. Although current models are numerically well defined and are able to resemble the tissue differentiation and adaptation processes, they are limited by (1) the fact that some of their input parameters are likely to be site-and species-dependent, and (2) their verification is done by data that may make the model results redundant. However, some theories do have predictive power despite the limitations of generalization. It seems to be a matter of time until new experiments and models appear with predictive power and where rigorous verification can be performed.


Abstract
It is well known that the mechanical environment affects biological tissues. The importance of theories and models that aim at explaining the role of the mechanical stimuli in process such as differentiation and adaptation of tissues is highlighted because if those theories can explain the tissue's response to mechanical loading and to its environment, it becomes possible to predict the consequences of mechanical stimuli on growth, adaptation and aging of tissues. This review aims to present an overview of the various theories and models on tissue differentiation and adaptation, and their mathematical implementation. Although current models are numerically well defined and are able to resemble the tissue differentiation and adaptation processes, they are limited by (1) the fact that some of their input parameters are likely to be site-and speciesdependent, and (2) their verification is done by data that may make the model results redundant. However, some theories do have predictive power despite the limitations of generalization. It seems to be a matter of time until new experiments and models appear with predictive power and where rigorous verification can be performed.

Introduction
It is well known that the mechanical environment affects biological tissues, in which mechanical loads activate or inhibit processes involving the genesis, adaptation, and aging of tissues [1]. This is especially important for tissues that have a primary biomechanical function such as musculoskeletal or cardiovascular tissues; however, biomechanical and mechanobiological factors appear to be critical for regulating cell behavior, and tissue maintenance and transformation in virtually all other tissues of the body [2].
The study of tissue differentiation and adaptation to their mechanical environment is challenging. These processes are complex, involving so many variables that physical experimentation is often either time consuming, expensive, or impossible [3]. For instance, forces applied to tissue might be of different type or nature: internal quasi-static forces caused by tissue growth; external forces imposed on the organism; and intermittent joint forces caused by muscle contractions [4]. As a result of the mechanical environment created by all those forces, time-dependent, spatially complex patterns of internal tissue stresses and strains are created in all tissues [4].
Since the 19 th century, several theories and models have been proposed to explain the role of mechanical stimuli on biological processes [5]- [13]. In the last decades, the increasing computational power has allowed the implementation and simulation of many of the proposed theories and models of tissue differentiation and adaptation. In mechanobiology, a field that describes the mechanisms by which mechanical loads regulate biological processes, computational models have been developed to complement lab experiments in order to propose and test rules that may govern the effects of mechanical loading on cells and tissues [14].
The importance of such theories and models is highlighted because if those theories can explain the tissues' response to mechanical loading, it becomes possible to predict the consequences of mechanical stimuli on growth, adaptation, and aging. Additionally, it is now clear that mechanobiological interactions between cells and their underlying substrate can critically influence cell behavior and the resulting tissue, even in tissues and organs without a predominant biomechanical role [2], [15]. Currently, models of tissue differentiation and adaptation gain importance due to the studies on functional tissue engineering, in which the effects of biomechanical factors on native and repair tissues, and the development and application of computational models of tissue growth and remodeling are among its priorities [15].
Several works have reviewed the processes of tissue differentiation and adaptation in musculoskeletal tissues from different perspectives [16]- [18]. However, to the knowledge of the author, the literature that reviews their mathematical fundamentals and implementation is scarce. For example [19], this review aims to present an overview of the various theories on tissue differentiation and adaptation, their mathematical implementation, and a discussion and possible future steps in this research area.

Classic Theories
In the late 19 th century, W. Roux (1981) suggested a relationship between mechanical stimuli and tissue differentiation [5]. His thesis, the "Theory of Functional Adaptation," was that cells are in constant competition among themselves for their own "functional stimulus" and, eventually, this competition determines the final tissue phenotype. The "functional" stimuli were defined as compression for the formation and maintenance of bone, tension for connective tissue, and a combination of both with shear for cartilage. Like Roux, F. Pauwels based his theory almost exclusively on clinical observations and simple experiments [11]. However, he introduced basic concepts of stress and strain to the description of the problem. In his theory, he identified the mechanical stimuli that are "felt" by the cells as elongation and hydrostatic pressure ( Figure 1). Pauwels proposed elongation as the stimulus for the formation of connective tissue and hydrostatic pressure as the stimulus for the formation of cartilage. In this theory, bone formation does not have a specific stimulus and only occurs by calcification of the cartilage when "the cell is protected from all direct mechanical stresses" and under the influence of hydrostatic pressure, see Figure 2. With the consideration of elongation and hydrostatic pressure as stimuli, Pauwels effectively described the tissue differentiation process in terms of mechanical invariants, i.e., magnitudes independent of the coordinate system.

Interfragmentary Strain Theory
Considering that strain tolerance is higher for granulation tissue, moderate for cartilage, and lower for bone tissue, Perren and Cordey [6] proposed that at the beginning of the fracture-healing process, when motion and strain are high, the first tissue that appears is granulation tissue since it is the only one able to withstand this mechanical environment. Once the granulation tissue stiffens the mechanical environment, the formation of cartilage can begin. Finally, with the presence of cartilage inside the gap, the mechanical environment becomes more suitable for bone generation. Perren and Cordey defined the concept of Interfragmentary Strain (IFS) for fractures as: where L is the length of the "fracture" gap and ΔL is the change in length caused by loading. Both terms are defined in the axial direction of the bone. The IFS is a measurement of stretch capacity that can be related to adequate environments for specific tissues.

Consideration of the Load History
Carter et al. were the first to explicitly interpret Pauwels' theory in terms of welldefined mechanical magnitudes and to include loading history in the analysis. They associated elongation with either octahedral shear strain or octahedral shear stress, and hydrostatic pressure with either hydrostatic stress or volumetric strain [21]. The hypothesis of Carter et al. involves these invariants with their load history, assuming that the variation in time of the mechanical loading is the triggering factor that stimulates tissue differentiation. Finally, Carter et al. proposed the use of an index that measures the load history of the suggested invariants. This index, initially considered only for osteogenesis, was called the Osteogenic Index (OI) [7]; however, a similar combination of hydrostatic stress and octahedral shear stress histories could be used for other tissues, see Figure 3.

The Fluid Flow as Stimulus
Considering that there is evidence showing the importance of fluid flow at the cellular level [22], Prendergast et al. described how fluid flow can amplify cellular deformation [9], [23] and thus proposed the use of the relative velocity between the phases as a stimulus for tissue differentiation. The biphasic theory developed by Mow [24] was used to determine the substrate strain as well as the relative velocity between the phases. Source: taken from Carter [21] Prendergast et al.'s theory affirms that the mechanical environment, in terms of the substrate strain and the relative fluid velocity, is the relevant factor in tissue differentiation. Finally, they hypothesized that changes in the mechanical environment are systematic and, consequently, that the mechanical stimuli could stimulate the replacement of one cell population by another. This means that the stimuli are not only the key factor for the generation of tissues, but are also important for the constant adaptation of tissues [9]. This theory is illustrated in Figure 4.

A Simple Quantitative Theory
Claes and Heigele [8] followed the main idea of Pauwels and Carter et al. in considering that the most important mechanical stimuli are strain and hydrostatic pressure. However, they included clinical information on ossification paths, which may guide the fracture-healing process. Their hypothesis was that "new bone formation in fracture healing occurs primarily along fronts of existing bone or calcified tissue and that the type of bone healing […] depends on the local strain and stress magnitudes". They used a finite analysis model to estimate the local magnitudes of strain and stress in callus tissues and along different ossification paths, which were obtained from histological images of a previous in vivo experiment.
The analysis performed by Claes and Heigele was a visual comparison of the typical locations of intramembranous bone formation and endochondral ossifications with the calculated strain and hydrostatic pressure. This analysis allowed them to state the quantitative relationships between mechanical stimuli and tissue differentiation shown in Figure 5. These relationships, determined by the mechanical stimuli, define specific areas in the figure that correspond to endochondral versus intramembranous ossification.

The Strain as Unique Mechanic Factor in Tissue Differentiation
In 2005, Duda et al. [25] presented a model of osteochondral defect healing that is based on the hypothesis that the strain of the tissue, expressed as the minimum principal strain, is a main stimulus for tissue differentiation. However, Duda et al. did not consider local stress as an additional stimulus, as Claes and Heigele did. Despite the simplicity of the model, it showed good qualitative and quantitative agreement with the histological findings.

Other Factors Different than Mechanical Loading
Biochemical signaling, cellular mechanisms, and local vascularization also play important roles in tissue differentiation and adaptation. Therefore, recent works have developed algorithms that emulate the effects of those factors, sometimes including them in the models of mechanically induced adaptation, in the search for a more comprehensive representation of the differentiation and adaptation processes.
Bailon-Plaza and Van der Meulen [26] proposed an early mathematical framework in which the effect of growth factors, as biochemical signals, was included in the process of bone healing. This model was later extended to include mechanically induced adaptation [27]. Isaksson et al. [10] and Moore et al. [28] also combined biological and mechanical factors in a mechanistic model of bone healing. They accounted for cell activities such as migration, proliferation, differentiation, and apoptosis. Quantification of the model constants, which mainly dictate the rates of various processes, was seen as a problem because they likely depend on the conditions of the subject (e.g., species, gender, clinical status, age, etc.).
Additionally, vascularization plays an important role in tissue generation and maintenance. To test how this factor may influence bone healing, Geris et al. [29] include the possibility of angiogenesis in the model previously developed by Bailon-Plaza and Van der Meulen [26]. This was done by describing the tensile Source: developed by Claes and Heigele [8] angiogenesis as a process involving endothelial cells, a vascular matrix, and a generic angiogenic growth factor produced by osteoblasts and hypertrophic chondrocytes.

Mathematical Implementation
The first implementations of the theories of tissue differentiation and adaptation used static mechanical models to determine the general stress state during a specific moment of the biological process. Currently, the models for the various theories have improved through the use of finite element analysis within an iterative algorithm that allows the study of the complete development of the process. Recent studies have commonly included dynamic variation of material properties, cell distribution, and mechanical loading. All the algorithms used have similar outlines (see Figure 6). First, there is an initial step in which the characteristics of the tissue are established (i.e., geometry and material properties). Second, the tissue is mechanically loaded and biophysical stimuli are determined (e.g., stresses, strains, flow velocities). Third, the phenotype for each point in the geometry is updated based on a numerical approximation of a chosen theory. Fourth, if the migration of pluripotent cells is considered, a new cell distribution within the tissue is obtained using, for instance, diffusion analysis. Fifth, the material properties are updated, and a new iteration, corresponding to a new time point, can be started. Convergence of the numerical procedures is often evaluated in terms of the mechanical properties; that is, if the change in mechanical properties is lower than a certain limit, the algorithm converges.
In a numerical description of the tissue differentiation and adaptation processes, it is possible to identify the following three different parts or (sub-) models: • A mechanical (sub-) model, which calculates the mechanical state of the tissue: stress/strain, flow velocities, and pressures. • A diffusion (sub-) model, which estimates the motion of cells and bio-signals due to concentration gradients. • A reaction (sub-) model, which predicts the change in the phenotype of the tissue cells as a result of the mechanical state of the tissue and the concentration of bio-factors. These are the numerical implementations of the previously described theories.
Brief descriptions of each model and their numerical implementations are presented below.

Mechanical Models
Determining the mechanical state of a tissue requires a proper description of the tissue as a material. Frequently, the tissues in the previously described theories are modeled as single-phase elastic and biphasic materials. The choice between these two material descriptions relies on the loading and boundary conditions as well as on the level of consideration of the model. In general terms, a single-phase elastic model is usually appropriate when strain rates are higher than 0.1 or 0.01 Hz, when the boundary conditions constrain free fluid flow from and to the tissue (low permeability), or when the consideration level is larger than the tissue level. Otherwise, a biphasic description should be used [30].

Single-Phase Elastic Model
Considering isotropy and linear elasticity (Hooke's law), the equations that relate stress and strain are as follows: with D e as the elastic stiffness matrix.
Finally, invariants such as the hydrostatic stress (σ d ) and octahedral shear stress (σ s ), parameters for the other (sub-) models (diffusion and reactive models), are calculated as follows:

Biphasic Approach
A more accurate material description is given by the biphasic theory or by poroelasticity. When incompressibility is assumed, the mathematical implementations and solutions of these techniques are equivalent [31], [32], and both approaches estimate the stress state of the phases as follows [33]: , where σ s is solid stress, σ f is fluid stress, e and ε are dilatational strain and total strain, respectively; p is apparent pressure, φ is the volume fraction, and λ and μ are the Lamé constants.

Diffusion Model [34]
Differentiation and adaptation of tissues are influenced by the local concentration of cells and bio-signals (e.g., growth factors). Motion of both, cells and biosignals, can be modeled as following a diffusion process. However, this diffusion process can have either of two approaches: the chemical pattern approach and the mechanochemical approach. On the one hand, in the chemical pattern approach, pattern formation and morphogenesis take place sequentially. Morphogenesis, i.e., formation of the structure, is essentially a slave process that is determined once the chemical pattern has been solved. On the other hand, in the mechanochemical approach, pattern formation and morphogenesis are considered to go on simultaneously as a single process, which may be closer to reality. The diffusion model that uses the mechanochemical approach consists of three governing equations: The conservation equation for the cell population density, the mechanical balance of the forces between the cells and the extracellular matrix (ECM), and the conservation law governing the ECM. Without minor details, the equations are written as follows: where the terms from left to right in the right side of the equation are related to cell convection, diffusion, haptotaxis, and mitosis, respectively; D 1 , D 2 , a 1 , a 2 , r, and N are positive parameters. The vector u is the displacement of ECM, and n and ρ are the number of cells and density of the ECM, respectively.

Cell-Matrix Mechanical Interaction Equation
where F is the external force acting on the matrix and σ is the stress tensor. However, the stress tensor has contributions from the ECM and the cells, such that: where σ ECM as a linear viscoelastic material, can be written in terms of the viscous and elastic components as: The subscript t denotes partial differentiation, I is the unit tensor, and μ 1 and μ 2 are the shear and bulk viscosities of the ECM. ε is the strain tensor, defined as , and θ is the dilation, On the other hand, the stress due to the cell traction, σ CELLS , is: where τ is the traction force generated by the cell, λ is a measure of how the force is reduced because of neighboring cells, and γ is a measure of the nonlocal long-range cell-ECM interactions. Finally, to complete the description of the terms in the cell matrix: where s > 0 is an elastic parameter characterizing the substrate attachments.

Matrix conservation equation
The conservation equation for the matrix material, ρ(r,t), is where matrix flux is considered due to convection and S(n,ρ,u) is the rate of secretion of matrix by the cells.

Reaction (Sub-) Model
The theories described in the first section can be written as numerical algorithms to be included, along with the mechanical and diffusion (sub-) models, in a model of tissue development and adaptation. Various authors highlighted above, among others, have implemented their theories using different mechanical stimuli, mechanical thresholds, or the methods of modeling the tissue. The mathematical implementations are described briefly in this section.

The Osteogenic Index
Carter et al.'s theory is qualitative, and no specific relationships between the stimuli have been proposed to define the regions shown in Figure 9. However, to achieve physiologic ossification patterns in their models, Carter  where a cyclic load condition has been discretized in n consecutive loading cases. For each point of the material one IO maxmin can be found using the maximum shear stress and the minimum value of hydrostatic stress during the entire cycle load [35]. The value of constant k must be previously calculated for each specific case (specific geometry, specific load conditions). This parameter is calculated by iteration until the plot of the OI matches the documented ossification patterns. Finally, matching plots of strain and hydrostatic stresses with histological images allows the determination of levels of pressure and strain for each phenotype. The various models of Carter et al. have used an elastic material model to represent the tissue, in a single-phase approach. This is held by the authors with the premise that the solution is sufficient to describe the ossification patterns, while other problems, such as the behavior of the cartilage, must be handled with biphasic theory [36].

Biphasic Implementation
When authors want to consider the relative flow as a stimulus for tissue differentiation, they must apply a biphasic theory. In the constitutive equations it is possible to relate the boundary conditions with the internal mechanical loads such as substrate strain, poro-pressure and/or relative velocity among phases, and the material properties of the tissue. Huiskes et al. [37] proposed some rules for tissue differentiation using the biphasic approach that were later refined by Lacroix et al. [38], and Kelly and Prendergast [39] as follows: γ + > for fibrous connective tissue, 1 3 v a b γ < + < for fibrocartilaginous tissue,  [40]. However, these values should be used cautiously because these constants are likely site-and species-dependent.

Discussion and Future Steps
The first theories of tissue differentiation and adaptation were based on simple experiments and clinical experiences. Later, similar hypotheses included mechanical variables in the description of the process. The implementations of mechanical variables and computational power have allowed not only the description, but also the prediction of events. Currently, models include other variables and parameters such as biosignaling and vascularization to increase their prediction capacity. Although the nature of these recently-included variables is distinct from the mechanical state of the tissue, they are influenced by the mechanical environment of the tissue and its function by being subject.
However, although it is accepted that mechanical function can influence the shape and structure of skeletal tissues, as well as the regulation of their development, it is still not totally clear how these processes are accomplished and how different inducing factors work together. Experimental verification of the different theories is out of the scope of this review. Nevertheless, verification is always an essential part of modeling, and the theories presented here have an important setback: their models often need input data that obviates the need for the models in the first place [41].
In general, tissue regulation includes two biological processes: an initial tissue differentiation and a subsequent tissue adaptation. The boundary between these processes is not well defined. To study both process using a single theory or model would require a common definition of change-inducing factors. However, the specific weights of those common change-inducing factors on the overall process are likely to change over time and based on tissue maturity, making it difficult to achieve a unified theory of tissue regulation.
The different quantitative relationships in the reaction models have been widely tested using computer models [29], [39]. However, conclusions from those works are always limited by the fact that many of the parameters have been defined based on specific experiments. This imposes a limitation because the value of such parameters is very likely to be site-and species-dependent [42], [43]. Moreover, the data needed to validate mechanoregulation theories often provide the same insight as the models, making the models less useful in certain cases of basic research [41]. However, some theories do have predictive power despite the limitations of generalization. It seems to be a matter of time until new experiments and models appear with predictive power and where rigorous verification can be performed [44].
Current theories and models are proper steps toward a final model that remains elusive. In general, numerical models force us to describe every detail of a theory in a quantitative way and, at the same time, allow us to make predictions instead of merely describing the process. However, the need for new experiments to evaluate the different tissue differentiation theories is evident. Testing of more and more complex research hypotheses is encouraged by the new and more powerful computational resources, but the original theories still require definitive validation or verification.