article: Mechanical stimulation of the healing tendon is thought to regulate scar anisotropy and strength and is relatively easy to modulate through physical therapy. However, in vivo studies of various loading protocols in animal models have produced mixed results. To integrate and better understand the available data, we developed a multiscale model of rat Achilles tendon healing that incorporates the effect of changes in the mechanical environment on fibroblast behavior, collagen deposition, and scar formation. We modified an OpenSim model of the rat right hindlimb to estimate physiologic strains in the lateral/medial gastrocnemius and soleus musculo-tendon units during loading and unloading conditions. We used the tendon strains as inputs to a thermodynamic model of stress fiber dynamics that predicts fibroblast alignment, and to determine local collagen synthesis rates according to a response curve derived from in vitro studies. We then used an agent-based model (ABM) of scar formation to integrate these cell-level responses and predict tissue-level collagen alignment and content. We compared our model predictions to experimental data from ten different studies. We found that a single set of cellular response curves can explain features of observed tendon healing across a wide array of reported experiments in rats–including the paradoxical finding that repairing transected tendon reverses the effect of loading on alignment–without fitting model parameters to any data from those experiments. The key to these successful predictions was simulating the specific loading and surgical protocols to predict tissue-level strains, which then guided cellular behaviors according to response curves based on in vitro experiments. Our model results provide a potential explanation for the highly variable responses to mechanical loading reported in the tendon healing literature and may be useful in guiding the design of future experiments and interventions.
  Many mechanically loaded tissues including skin, tendon, ligament, and heart respond to injury by forming a collagen-rich scar. Normally, tendons are comprised of highly aligned collagen fibers that help transmit forces between muscles and bones throughout the body and bear high loads. The Achilles tendon in particular can be exposed to loads up to 70 MPa, compared to 30 MPa in most other tendons [1]. These high loads often lead to injury, with Achilles tendon ruptures accounting for up to 45% of all tendon ruptures [2] and afflicting up to 2. 5 million annually [3,4]. Many who suffer from an Achilles tendon rupture never regain complete function, especially because healing tendons form scar with reduced collagen fiber organization and stiffness compared to uninjured tendons [5]. While there seems to be general agreement that mechanical stimulation of the healing tendon, such as during physical therapy, influences scar mechanical properties [6], current treatments for patients with an Achilles tendon rupture have produced variable results [1,3, 7–9]. To better understand the impact of loading, rat animal models of Achilles tendon rupture have been utilized so that scar tissue can be excised to determine quantitative biomechanical scar properties. Unfortunately, these studies have also led to a wide variety of results, with mechanical loading sometimes appearing to increase, but at other times appearing to decrease, tendon properties such as stiffness or rupture strength. One of the potential explanations for this variability could be due to differences in mechanics during healing, which could be altered through unloading (e. g. cast immobilization) or loading (e. g. free cage activity) of the tendon [10–14] as well as through the choice to repair the tendon with suture versus allowing natural healing without repair [12,15]. In this work, we developed a multiscale computational model of the healing rat Achilles tendon that integrates information about how local mechanics influences cellular alignment and collagen remodeling to predict the effects of various repair and loading protocols on tendon structure. We found that the multiscale model predicted the major observed trends in the evolution of tissue-level scar properties across a wide variety of published rat studies. Furthermore, the model simulations identified a potential mechanism underlying the apparently paradoxical finding that mechanical loading enhances collagen alignment in unrepaired Achilles tendons yet decreases it in repaired tendons.
  We searched the literature for studies of transected Achilles tendons in rats, the most common animal model used to mimic an Achilles tendon rupture. These studies used a variety of time courses and mechanical loading protocols to treat the rats during the wound healing process. First, we focused on studies that tested the effects of natural healing of unrepaired, transected tendons and imposed unloading by either botulinum toxin (Botox) injection into the gastrocnemius muscle or tail suspension [10,24–26] or loading by allowing the rat to freely walk around its cage [10,13,27,28] (Fig 2A). Because no studies quantified both collagen alignment and total collagen mass, we then narrowed our search to include only studies that reported at least one of two quantitative measures that could act as surrogates for model comparisons. Based on the strong reported correlation between measured collagen alignment and intrinsic material properties of the tissue [6,14], we compared measured values of Young’s modulus (Ey in MPa) to the levels of collagen alignment predicted by our model. Similarly, based on previous studies showing that collagen concentration and tendon CSA rise in parallel [24,29,30], as well as the fact that tissue mass increases with tissue volume, we compared measured values of tendon cross-sectional area (CSA in mm2) to the total collagen content predicted by the model. Despite considerable variability among the reported values of these metrics, two clear trends were apparent in the data. First, Ey started near 0 MPa at 3 days and rose to about 20–40 MPa at 14 days in both unloaded and loaded conditions (Fig 2B and 2C). Second, CSA values in unloaded conditions remained at around 5 mm2 (Fig 2E), while loading increased CSA over time to around 10–20 mm2 at day 14 (Fig 2F). Next, we used the model of the rat hindlimb implemented in OpenSim [16] to simulate these unloading and loading conditions and estimate strains in the healing region. We simulated unloading by fixing all joint ankles in plantar flexion, with minimum muscle activation and a 7mm gap distance between the tendon stumps [10,25,31], resulting in a predicted constant cell strain of E11 = 0. 002 in the primary loading direction (Fig 3A). We simulated loading by prescribing joint angles and muscle activation corresponding to the rat gait cycle. Predicted strains oscillated between. 009 and. 043, reflecting the swing/stance phases of the gait cycle and yielding a cyclic strain amplitude of 0. 034 and a mean strain of 0. 030 (Fig 3A). Using our published model of stress fiber dynamics [18], both strain states produced cell alignment along the loading axis, with slightly lower predicted alignment for the loaded case (Fig 3B). Using the collagen synthesis curve we fitted to published data [32–35], mean strain from the loaded case was associated with about 2x greater collagen production than the mean strain from the unloaded case (Fig 3C). Integrating these two cellular behaviors in the ABM component, the multiscale model predicted that the specific loading protocols we simulated should produce little difference in collagen alignment (Fig 2D) but a substantial difference in collagen content (Fig 2G). Trends in predicted alignment agreed with literature reports showing similar Ey for both groups at all time points (Fig 2B and 2C). Furthermore, differences in model-predicted collagen accumulation in the two loading states qualitatively matched reported differences of the tendon CSA in loaded (Fig 2F) compared to unloaded (Fig 2E) conditions. In our next set of simulations, we explored how surgically repairing rat Achilles tendons would alter our predictions and the response to loading during healing (Fig 4A). We selected a set of studies that subjected experimental groups to either unloading with cast immobilization or loading with an exercise protocol, defined as 60 min/day treadmill exercise + free cage activity for the rest of the time, compared surgically repaired and unrepaired groups, and reported Young’s modulus or cross-sectional area (Fig 4B and 4D) [12,13,36]. For the unrepaired groups in this second set of studies, trends were similar to those shown in Fig 2, with loading enhancing CSA but not Ey. Similar to the unrepaired groups, loading of suture-repaired tendons increased tendon CSA compared to unloading (Fig 4D). However, loading in the suture-repaired groups surprisingly reduced Young’s modulus below the values observed in any other group by the 42-day final time point (Fig 4B). Assuming that surgical repair eliminated the gap distance between the tendon stumps, our hindlimb simulations in OpenSim predicted higher strains for repaired conditions compared to unrepaired, with repaired-unloaded tendons experiencing a static strain of E11 = 0. 012 and repaired-loaded tendons oscillating between 0. 011 and 0. 075 (cyclic strain of 0. 064, mean strain of 0. 049; Fig 5A). The much larger cyclic strains in the repaired-loaded group induced stress fiber disassembly along the loading axis in the cell alignment model, resulting in cells that were nearly randomly oriented (Fig 5B); the differences in predicted stress fiber orientation distributions for these four cases are shown in Fig 6. On the other hand, the large mean strains in the repaired-loaded group led to a higher rate of collagen synthesis in this group compared to the other conditions simulated (Fig 5C). Integrating these predictions using the ABM component, overall our multiscale model correctly matched the apparently paradoxical reports that repaired, loaded tendons have the lowest collagen alignment of any of these four conditions (Fig 4B and 4C) despite having the highest collagen content (Fig 4D and 4E). In other words, a single set of cellular response curves can explain features of observed tendon healing across a wide array of reported experiments in rats without fitting model parameters to any data from those experiments. Rather, the key to the predictions in our multiscale model is simulating the specific loading and surgical protocols to predict tissue-level strains, which then guide cellular behaviors according to response curves based on in vitro experiments. Many of the studies we simulated here loaded healing tendons through unrestricted cage activity or through daily exercise added to normal cage activity. In these protocols, tendons are cyclically loaded in short bursts as the animals move about their cages, interspersed with short rest periods when they are standing still and longer rest periods when they sleep. We chose to model this situation by imposing cyclic loading at 1 Hz with a 1-hour on, 1-hour off protocol for 12 hours, followed by 12 hours of rest, and then repeating. To understand the effect of this choice on our results, we used the cell alignment model to simulate the effect of stretching with different protocols that all produced the same time-averaged frequency of 0. 5Hz (Fig 7A). We compared stretching for repetitions of 6 hours at 1Hz followed by 6 hours of rest (0Hz), repetitions of 1 hour at 1Hz followed by 1 hour at rest, and continuous stretching at 0. 5Hz. We simulated each repetition until steady state was reached, defined as the point when the difference in alignment at the end of two consecutive repetitions was less than 0. 01 (Fig 7B). We found that the different stretching protocols produced different relationships between strain amplitude and predicted steady-state cell alignment, defined as the average order parameter over the last simulated cycle. At low strain amplitudes, all the protocols produced similar results, while at peak strains between 0. 06 and 0. 10 the three protocols resulted in predictions ranging from fairly strong alignment parallel to stretch (continuous stretching) to a modest degree of perpendicular alignment (1h on / 1h off and 6h on / 6h off). These differences arise from two features of the stress fiber model. First, stress fiber disassembly is triggered by high negative strain rates in the model, so much higher strains are required to influence alignment when loading is imposed at a lower frequency (Fig 7A). Second, because disassembly of stress fibers occurs on a much faster time scale than assembly, shorter durations of cyclic stretching can drive alignment down quickly, while much longer rest periods are needed to recover from each loading cycle (Fig 7B). Our simulations of intermittent loading protocols raise the possibility that relatively minor differences in loading could have significant implications for fibroblast and collagen alignment. This could be an intriguing explanation for the surprising degree of variability we found in the literature among studies that employed apparently identical experimental protocols and outcome metrics (see Fig 2B, 2C and 2F), or even within groups in individual studies (see error bars in Fig 4B and 4D). To date, most published experiments have not tracked movement of rats with the detail required to simulate more realistic or even animal-specific loading protocols. Sams-Dodd observed that healthy rats traveled around 5700 cm within a 10 minute observation period in one of the only studies we could find that attempted to track travel distances [37]. Our simulations suggest that employing continuous movement tracking in future studies might provide additional insight into whether individual variability in activity and loading can explain some of the observed variability in tendon healing. The most serious limitation of the modeling studies reported here is that we were only able to validate them through qualitative comparisons of model-predicted trends to experimentally measured surrogates. The most novel prediction of our model was the degree of collagen fiber alignment. Experimentally, collagen alignment can be directly measured from ultrasound [4,12] or polarized microscopy [38]. However, these measurements were reported in so few studies that we were forced to use a more commonly reported surrogate, the Young’s modulus (Ey), which has been shown to correlate with the degree of alignment. It is conceivable that other factors such as collagen density that might differ between the cases simulated here might have influenced Ey. Future studies may be able to draw on relationships between alignment and modulus such as those reported by Lake et al. [39] and Li et al. [40] to quantitatively estimate Ey. The comparison of total collagen content to tendon cross-sectional area (CSA) could also have limitations, since experimentally the rate of collagen production influences both collagen density (as assessed by biochemical assays or picrosirius red staining) and total CSA, but few studies provide data on both simultaneously. In addition, our collagen content predictions are in arbitrary units, since the studies we used to determine the effect of stretch on collagen synthesis reported relative changes rather than absolute synthesis rates. The simplifications we made in the rat hindlimb model could also introduce some errors in the tendon-level strains we calculated. First, we included only the musculo-tendon units that directly comprise the Achilles tendon. Second, due to a lack of information in the literature about hip and knee motion during walking in the setting of rat Achilles tendon injury, we combined healthy hip and knee with injured ankle motion [41,42]. However, we expect errors introduced by this choice to be small, since the gastrocnemius and soleus muscles do not cross the hip, the soleus does not cross the knee, and the moment arms around the ankle are larger than at the knee for the gastrocnemius. Furthermore, as long as simulated strains from repaired tendons are higher than unrepaired and strains from loading are higher than from unloading, the overall trends predicted from our model should be robust to small changes in the exact magnitude of the predicted tendon strains. Our agent-based model assumes that cells deposit collagen aligned with their own axis, and this assumption was critical to translating cell alignment predictions from the cytoskeletal model into tissue-level predictions of collagen structure. While the exact mechanisms by which fibroblasts deposit and orient collagen in vivo are still being debated, the general idea that collagen ends up locally aligned with the fibroblasts that deposit it remains strongly supported in the literature [43–45]. Furthermore, we have previously shown that agent-based models incorporating this same assumption correctly predict a range of scar structures observed following myocardial infarction in various animal models under different mechanical conditions [23]. We also considered and simulated several other alternative methods of determining collagen alignment in the model but found that none could predict all the observed trends apparent in the data reviewed here. For instance, several experiments have theorized that surrounding collagen fibers could “structurally constrain” the formation of new fibers in vivo [46–48]. Others have demonstrated strain-dependent modulation of collagen degradation that could influence overall alignment under uniaxial loading by selectively degrading fibers with certain orientations faster than others [49–51]. While all these effects could be present within the actual tendon, we found that strain-dependent cell alignment, deposition of collagen aligned with the cells, and strain-dependent changes in collagen synthesis rate were sufficient for capturing the major trends in the data as outlined above. In this study, we used multiscale modeling to integrate information from the literature on fibroblast responses to stretch (alignment and collagen synthesis), scar formation following injury (collagen deposition and other features of the agent-based model), and musculoskeletal mechanics (rat hindlimb model implemented in OpenSim) to interpret apparently conflicting data from a range of experimental studies. We found that our computational model could reproduce several key features of observed tendon healing across a wide array of reported experiments in rats–including the paradoxical finding that repairing transected tendon reverses the effect of loading on alignment–without fitting model parameters to any data from those experiments. Rather, the key to the predictions in our multiscale model was simulating the specific loading and surgical protocols to predict tissue-level strains, which then guided cellular behaviors according to response curves based on in vitro experiments. These results suggest that the apparently conflicting data in the studies we reviewed may in fact reflect consistent biologic responses to local strains in the healing tendon, providing a new conceptual framework for interpreting existing data and devising potential therapies for Achilles tendon rupture.
  We searched PubMed and Google Scholar for all papers that included the keywords “rat Achilles tendon rupture injury” in the title or abstract and used a scalpel to perform a full transection of both the Achilles tendon and the plantaris tendon, a tendon parallel to the Achilles that is proportionally larger in rats than in humans and can act as an “internal splint” [52]. From these, we next identified studies that specifically compared the effects of unloading from intramuscular injection of Botox, tail suspension, or cast immobilization against loading due to free cage activity or treadmill exercise. From these 30 identified papers, we then selected those that quantitatively measured the cross sectional area and/or the Young’s modulus of the healing scar, leaving us with 10 studies that met our search criteria [10–14,24–26,36,53]. We used cross sectional area (CSA) of the tendon as a surrogate measure for collagen production, based on previous studies showing that collagen concentration and tendon CSA rise in parallel [24,29,30], as well as the fact that tissue mass increases with tissue volume. Similarly, we used the Young’s modulus as a surrogate measure for collagen alignment based on the strong reported correlation between measured alignment and intrinsic material properties of the tissue [6,14], and the fact that Young’s modulus was commonly measured in the studies we found while collagen alignment was not. We adapted a previously published rat right hindlimb model to conduct our simulations (Fig 1A) [16] (https: //simtk. org/projects/rat_hlimb_model). Briefly, the model used anatomically accurate representations of the bones (spine, hip, femur, tibia, and foot), joints, and musculo-tendon units (each represented by one line segment from its origin to insertion) of the rat right hindlimb. Musculo-tendon units were represented as linear elements consisting of a muscle segment in series with a passive tendon segment. The muscle segment consisted of its own passive element in parallel with an active contractile element that generated force depending on a Hill-type model. The mechanical properties of muscle fiber and tendon were defined using fiber force-length, fiber force-velocity, and tendon force-strain curves determined by Millard et al. [54]. Muscle segments could be prescribed activation levels varying from 0 (rest) to 1 (full activation). Since muscles not attached to the Achilles tendon were irrelevant for this simulation, we simplified the model by only including the lateral and medial gastrocnemius and soleus muscles, the three musculo-tendon units that comprise the Achilles tendon. To predict Achilles tendon strains for the various cases simulated here, we prescribed experimentally measured joint angle profiles, muscle activation curves, and tendon mechanical properties as inputs, and obtained strain vs. time curves from forward dynamic simulations using custom MATLAB routines employed by two of the authors in a previous publication [55]. We used a computational model of stress fiber remodeling published previously by our group to estimate cell alignment behavior [18] (Fig 1B). The model represents the thermodynamics of stress fiber (SF) assembly and disassembly, capturing features such as the ability of tension to promote assembly by altering the free energy of bound actin subunits. Two specific features of the model are important for the predictions shown in this manuscript. First, on the time scale of individual loading cycles, large negative strain rates reduce stress fiber tension through the force-velocity behavior of myosin, promoting SF disassembly in the direction of stretch and net SF (and cell) orientation perpendicular to that stretch. On a longer time scale, the model assumes that cells can remodel the extracellular matrix and/or their attachments to the ECM to attain an equilibrium strain state that minimizes the sum of the energies associated with elastic stretch and the chemical potential of the bound and unbound actin subunits. This aspect of the model drives the response to mean boundary conditions, upon which cyclic responses are then superimposed. Here, we used the average and spread of predicted SF distributions as surrogate measures of cell alignment, the same approach we used in the original model validation against in vitro data. For a more complete description of the model and details on its validation, please see Chen et al. 2018 [18]. We adapted an agent based model (ABM) originally published by Rouillard and Holmes [23,62] for infarct healing to integrate cellular-level responses and predict the evolving tendon scar structure (Fig 1D). Table 1 lists all model parameters altered for these simulations, while Fig 9 shows a flowchart of a cell’s decision tree within the model. Rouillard and Holmes modeled fibroblasts as circular discs free to move in a square, two-dimensional space divided into 10-micron-square patches. Each patch contained information about the local collagen alignment and density. Fibroblasts could migrate, proliferate, undergo apoptosis, and remodel collagen. Fibroblast orientation guided fibroblast migration direction and deposition of collagen, and existing collagen fibers were degraded at a rate proportional to their local concentration. A local chemokine concentration gradient with a high concentration of chemokines within the wound area and a low chemokine concentration in the healthy tissue drove cell migration into the wound (Fig 10A). We made several modifications to this model to adapt it for Achilles tendon healing. We simulated healing of a rectangular wound area after complete transection or rupture (Fig 10B). The initial wound size, matrix structure, and collagen content within the wound depended on the healing condition (Table 1). An unrepaired transected tendon contained a low amount (0. 1%) of randomly aligned collagen fibers, mimicking the randomly aligned provisional matrix in the wound area, while a suture-repaired tendon contained some aligned collagen (0. 9%, alignment order parameter of 0. 4), mimicking aligned collagen fibers from the healthy tendon stumps. Fibroblasts migrated into the wound space (Fig 10B, 10C and 10D) from the two opposing sides adjacent to healthy tissue. Rouillard determined cell alignment from a series of phenomenologic equations that represented the alignment response to stretch, contact guidance from surrounding collagen, and chemokine gradients, as well as their integration. Here, we replaced the original phenomenologic relationship governing stretch-induced alignment with the cell alignment predicted from the stress fiber model, and employed a more recent equation for integrating across the alignment cues published by Richardson et al. [63]. Based on previous studies showing that cells upregulate collagen production after exposure to both static [64] and cyclic stretch [32–35], cells synthesized collagen according to their mean strains, with higher strains corresponding to higher collagen synthesis rates. The mean strain was calculated by taking an average over one gait cycle period (1 second). In unloaded cases, this mean strain matched the static strain values (Fig 5). We determined the collagen synthesis amounts according to a sigmoidal curve fitted to data from four independent experiments in the literature (Fig 1C) [32–35]: CollagenSynthesis (εm) =1. 31+exp (−150 (εm2−0. 02) ) +0. 6, (2) where εm is the mean strain felt by the cell. Cells deposited collagen aligned to their major axis of alignment. While the exact mechanisms by which fibroblasts deposit and orient collagen in vivo are still being debated, the general idea that collagen ends up locally aligned with the fibroblasts that deposit it remains strongly supported in the literature [43–45]. At each time point for which collagen content and orientation are reported, the local collagen content and orientation histograms from each collagen patch were averaged to determine a single area fraction, mean angle, and order parameter for the entire scar region.
summary: Tendons and ligaments transmit force between muscles and bones throughout the body and are comprised of highly aligned collagen fibers that help bear high loads. The Achilles tendon is exposed to exceptionally high loads and is prone to rupture. When damaged Achilles tendons heal, they typically have reduced strength and stiffness, and while most believe that appropriate physical therapy can help improve these mechanical properties, both clinical and animal studies of mechanical loading following injury have produced highly variable and somewhat disappointing results. To help better understand the effects of mechanical loading on tendon healing and potentially guide future therapies, we developed a computational model of rat Achilles tendon healing and showed that we could predict the main effects of different mechanical loading and surgical repair conditions reported across a wide range of published studies. Our model offers potential explanations for some surprising findings of prior studies and for the high variability observed in those studies and may prove useful in designing future therapies or experiments to test new therapies.
keywords: dermatology
  medicine and health sciences
  surgical repair
  vertebrates
  animals
  fibroblasts
  surgical and invasive medical procedures
  collagens
  simulation and modeling
  animal models
  mammals
  physiological processes
  model organisms
  systems science
  mathematics
  connective tissue cells
  experimental organism systems
  research and analysis methods
  computer and information sciences
  scars
  tissue repair
  animal cells
  proteins
  animal studies
  agent-based modeling
  connective tissue
  biological tissue
  biochemistry
  rodents
  eukaryota
  anatomy
  cell biology
  tendons
  physiology
  biology and life sciences
  cellular types
  physical sciences
  amniotes
  organisms
  rats
title: Multiscale computational model of Achilles tendon wound healing: Untangling the effects of repair and loading