A Coupled Mechanobiological Model of Muscle Regeneration In Cerebral Palsy

Introduction

Cerebral palsy (CP) is a neuromusculoskeletal disorder arising from a static neural lesion but leading to musculoskeletal and gait impairments that can give rise to long-term degradation of musculature (Fridén and Lieber, 2003; Smith et al., 2011; Larkin-Kaiser et al., 2019). CP is the most common cause of physical disability in children and manifests as spasticity, contractures, poor control of muscles, and altered reflexes and posture (Graham et al., 2016). Depending on disease severity, muscles in individuals with CP are often smaller and weaker than typically developing (TD) counterparts; additionally, muscle size and strength decline over time (Elder et al., 2003; Handsfield et al., 2016; von Walden et al., 2017; Sahrmann et al., 2019). The macroscale changes to CP muscle are well-known; however, cellular level studies of muscle regeneration are only commonly performed for TD muscle, as opposed to CP muscle. Of the few recent studies performed on CP muscle, cellular level differences between CP and TD muscle have been found and include increased collagen deposition in the extracellular matrix (ECM) (Booth et al., 2001; Fridén and Lieber, 2003; De Bruin et al., 2014), decreased number of muscle stem cells (Smith et al., 2013), decreased stem cell activity (Domenighetti et al., 2018), and an increase in pro-inflammatory gene expression compared to TD muscle (Von Walden et al., 2018).

Skeletal muscle is a post-mitotic tissue capable of repair and regeneration. Muscle regeneration and repair may be triggered by different cues including trauma, muscle use and strain, and chronic degenerative diseases, which over time lead to tissue adaptation (Artrong et al., 1991; Järvinen et al., 2008; Chen et al., 2020). Typically, eccentric exercise stimulation of muscles attenuates age-related muscle loss and promotes myofiber hypertrophy (Chen et al., 2020). Stimuli such as eccentric lengthening exercises cause mechanical strains in the muscle that damage cell membranes and lead to a cascade of chemical signals and cellular responses. Following damage, the muscle fiber environment undergoes a tightly regulated adaptive repair process which is often categorized according to a series of four phases of regeneration: 1) damage in the form of membrane rupture, 2) acute inflammatory response from macrophages and neutrophils, which involves breakdown and clearance of necrotic tissue, 3) regeneration orchestrated by activation, proliferation, differentiation of myogenic precursor cells and fusion of myoblasts to the debrided region of the myofiber, and 4) repair and remodeling of the ECM by fibroblasts (Partridge, 2002; Mourkioti and Rosenthal, 2005; Chargé and Rudnicki, 2009; Novak et al., 2014; Laumonier and Menetrey, 2016). During the first step, the breakdown and necrosis of myofibers is triggered via the disruption of the sarcolemma and subsequent increase in permeability. Creatine kinase is released into the serum and is a common marker of post-mechanical stress or muscle degeneration. Loss of calcium homeostasis and calcium influx, due to damage of the sarcolemma, then drives tissue necrosis. The result of injury is focal or total autolysis of fibers (Chargé and Rudnicki, 2009). Step two is marked by degeneration that occurs within hours of damage with the activation of myeloid and secretory cells that predominantly release cytokines. Neutrophils invade injured muscle within 1–6 h and reach peak concentration 6–24 h post-injury. At 48 h post-injury, macrophages become the predominant inflammatory cell types at the injury site. Macrophages phagocytose cellular debris and activate myogenic cells, ready for the regeneration process (Novak et al., 2014). Their importance in skeletal muscle regeneration is due to their phagocytic and antigen-presenting roles (Tidball and Villalta, 2010). Arnold et al., 2007 postulate that phagocytosis of muscle cell debris induces a switch of pro-inflammatory macrophages toward an anti-inflammatory phenotype, releasing TGF-β. This also suggests that inflammatory macrophages stimulate myogenic proliferation while anti-inflammatory macrophages exhibit differentiating activity.

The damaged environment activates quiescent satellite cells (SCs) and fibroblasts that remodel the affected tissue. SCs are a mononucleated, progenitor cell population pivotal to physiological muscle repair and regeneration. In an uninjured state, SCs sit between the plasma membrane of the muscle fiber and the basal lamina. SC content of muscle differs between age groups and activity levels. Increased SC density is observed at the neuromuscular junction and adjacent to capillaries. This suggests the muscle environment created by and surrounding these structures may attract SCs or regulate the distribution of the SC pool. The regulation of SC density is also observed down to the single myofiber level and during regeneration, activation of SCs is not restricted to the site of injury on the myofibre. Mechanical stimulation through endurance and resistance exercise can also accelerate the turnover of ECM components in skeletal muscle. The ECM is primarily composed of collagens, laminins, fibronectin and proteoglycans. Fibroblasts synthesize and assemble most of the ECM in skeletal muscle, while other components are responsible for degradation of the ECM (Lu et al., 2011). Collagen synthesis is increased post-exercise and transcriptional analysis shows increased encoding of collagen types I, III and IV post-endurance training (Grzelkowska-Kowalczyk, 2016). A loss of balance in terms of fibroblast secretion of ECM and clearance of collagen may result in perturbed muscle regeneration and fibrosis.

The muscle regeneration process builds new healthy muscle under optimal conditions and can be organized and considered as a four step process (Chargé and Rudnicki, 2009). It should be noted that perturbation of this process may lead to fatty infiltration and fibrotic tissue (Joe et al., 2010; Uezumi et al., 2010; Wang et al., 2015). While regeneration in skeletal muscle occurs at the cellular level, degeneration of CP muscle over time leads to both cellular and observable macroscale changes (Graham et al., 2016). In light of this, it bears considering whether the process of muscle regeneration may lead to degeneration in CP primarily as a result of changes to the cellular environment. The regeneration process is complex however, and exploration of this problem requires understanding interactions at multiple scales among multiple cellular agents.

Agent-based modeling (ABM) is a technique capable of probing complex adaptive systems from the bottom-up. In ABM, autonomous agents are situated in an environment with changeable relationships. ABM is well-suited for biology as bottom-up modeling enables cells to act and react to one another and to local stimuli without an a priori macroscale outcome. This is achieved through its representation of multiple levels of biological organization, capturing of intracellular dynamics between large populations of cells, and its ability to integrate cell-signaling events (An et al., 2009). Furthermore, both cellular and non-cellular components of an agent-based model can be programmed to perform biologically relevant behaviors such as proliferation, apoptosis and migration (Gorochowski, 2016). The macroscale behavior observed is then either a directed or an emergent effect of the local cellular actions and interactions; thus, ABM is a useful tool for capturing the complex, nonlinear, and multiscale nature of physiology (North et al., 2013; Wilensky and Rand, 2015), and is a promising approach to model the process of muscle regeneration. The ABM approach has previously been applied to studies of tumor formation, cardiac modeling, vascular remodeling, bone tissue regeneration, wound healing, signaling, and metabolic processes (Bailey et al., 2009; Flegg et al., 2015; Borgiani et al., 2017).

ABM has been used previously to explore muscle-related pathologies such as Duchenne muscular dystrophy (Virgilio et al., 2018) and disuse atrophy (Martin et al., 2015). Here, we develop an ABM that simulates typical muscle regeneration based on physiological properties and rules derived from literature. To simulate a muscle’s response to mechanical stimulus, we built a 3D finite element (FE) model of a muscle fiber bundle that underwent active eccentric contraction. The resultant areas of high strain and thus mechanical damage in the FE model serve as cues for remodeling in the agent-based model. This coupled model links mechanical stimulus and damage to its physiological response in skeletal muscle. We use this model to simulate the mechanobiological feedback loop of muscle repair over 3 months and investigate chronic regeneration and degeneration in TD and CP muscles. The purpose of this study was to simulate active eccentric contraction of muscle to obtain local strain data, then to seed the highest localized strain points into the agent-based model of muscle regeneration and evaluate the sensitivity of fibril recovery post-injury after altering SC concentration.

The pathological differences in muscle geometry and deformation that may lead to changes in the tissue microstructure can be observed by coupling ABM and FE modeling. Firstly, the FE analysis is used to determine areas of high strain and to localize tissue injury. These areas direct cell migration in the agent-based model as it is used to simulate a cycle of repair following injury. Secondly, following the repair time course, the resultant agent-based model geometry and cell counts can inform the reconstruction of the FE model geometry, thereby providing progressive morphological data. Lastly, as the two models provide repeated feedback, the effect of pathological muscle morphology on chronic injury and regeneration can be observed and compared to typical muscle regeneration and muscle morphology. By coupling FEM with ABM in an iterative fashion, the biological processes of tissue adaptation can be explored over time. This coupled model investigates how impairment in the muscle regeneration process can influence pathological muscle degeneration in CP. There is a need to better understand this process in individuals with CP over time.

Methods

This study coupled an agent-based model and FE model of a human skeletal muscle fiber bundle to simulate muscle regeneration in response to eccentric contraction in CP and TD muscle. Initial geometry for all models was built from a single human muscle fiber cross-section obtained from the literature (Mackey et al., 2017), where this group undertook antibody staining of muscle biopsies from four young healthy male human participants. The vastus lateralis is a frequently used muscle in biopsy studies since it is a large, easily accessed muscle (Baczynska et al., 2016). In this case, we used a histological image from human vastus lateralis as a representative of fiber bundle architecture for human lower limb muscle. The histological image of muscle fiber bundle cross-section was segmented to distinguish 20 muscle fibers from the extracellular matrix (ECM) (Figure 1) using the Statistics and Machine Learning Toolbox™ in MATLAB (The MathWorks, Inc., Natick, MA). The 20 muscle fibers form a fiber bundle with each fiber surrounded by ECM and allows the study of cell-cell and cell-ECM interactions. Briefly, color-based K-means clustering was used to distinguish ECM from fibers, labeling every pixel in the image with a cluster index (1 or 0). The following sections discuss the agent-based model and FE model construction separately.

FIGURE 1. Geometry for both the agent-based models and FE models were generated from a single histological section. The coordinates for each pixel were extracted in MATLAB and used to recreate (A) pixels on ABM grid, and (B) to extrude muscle fibers for the FE models.

Agent-Based Modeling

The agent-based model was implemented in Repast Simphony (North et al., 2013), a Java-based modeling platform. Pixel values of the initial geometry were loaded onto a grid at the corresponding coordinate points in Repast Simphony (Figure 1A). The agent-based model contains 20 muscle fibers and represents a cross-sectional slice thickness of 50 µm. The ABM rules were developed based on literature descriptions of physiological interactions (see subsections below) (Figure 2). The model was built to simulate the progression of events during muscle regeneration known to take place over 4 weeks following injury (Laumonier and Menetrey, 2016).

FIGURE 2. Overview of agent-based model steps. FE model resultant strains are mapped from each element center coordinate to the corresponding pixel coordinates of fibrils in the agent-based model. Fibers, SCs, neutrophils, macrophages, fibroblasts and ECM components work to regenerate damaged tissue.

Cell populations modeled as agents in the agent-based model include muscle fibers, macrophages and neutrophils, SCs, fibroblasts, and ECM components. Extracellular guidance cues influenced the behavior of agents. These extracellular guidance cues comprised the cytokines and growth factors IGF-1, TGF-ß, HGF, IL-6 and TNF-α (Eqs. 1–5). Cytokine levels were based on cell-type specific secretions or generalized cell secretions. The agent-based model simulated muscle regeneration over 28 days following 10% damage to muscle fibrils. Initial agent-based model geometry consisted of 11719 2D grid elements. The model simulated 28 days of regeneration with each tick representing 1 hour of cell activity. Each pixel represents 6.45 µm² of tissue. Each simulation was run 150 times, unless otherwise stated, and simulation results are shown as mean ± standard deviation (SD).

Two ABM environments were developed to represent the muscle milieu in (1) a healthy TD skeletal muscle and (2) a CP skeletal muscle. Typical SC concentration in young adult males and children is reported as ∼0.10 SC per fiber for 10 µm slice thickness (Verdijk et al., 2014; Snijders et al., 2015; Mackey et al., 2017). A decrease in SC concentration ∼60% has been found in muscle biopsies from children with CP compared to typical adolescents (Smith et al., 2013). Therefore, SC density for 10 µm slice thickness was 0.10 SC per fiber for the healthy skeletal muscle and 0.04 SC per fiber for the CP muscle, representing a 60% decrease. For 50 µm slice thickness, the equation for SC counts was $SC density for 10\mu m\times number of fibers\times \frac{slice thickness}{10}$. Initial SC counts were set to 0.1 × 20 × 5 = 10 for the TD muscle and 0.04 × 20 × 5 = 4 for CP muscle. We additionally conducted a sensitivity analysis to assess the impact of variation in SC count on the resulting fibril count over time.

Fibers and Extracellular Matrix

The creation of fiber and ECM geometry has been described above. Each agent was initialized to a coordinate point on a 2D grid based on the segmentation results. Damage was assigned to fiber and ECM regions in the agent-based model based on regions of high strain in the FE model (described in detail below). Following eccentric contractions, regions of strain above a certain threshold in muscle indicate the locations of injury (Best et al., 1995; Garrett, 1996). In the agent-based model, fibril agents with the highest 10% of strain values from FE modeling were regarded as damaged in the agent-based model. Neighbor fiber and ECM agents were also set to damaged, where neighbors were defined by searching the von Neumann neighborhoods of damaged fiber agents. This extended area of damage simulated necrotic tissue (Laumonier and Menetrey, 2016). Damaged ECM had collagen density set to a lower value and signaled for fibroblasts to initiate the repair process (Schoenrock et al., 2018). When average fibril count per fiber was repaired to pre-injury levels, fiber borders were re-formed and additional fibrils could be added, consistent with the “ghost fiber” phenomenon previously observed (Webster et al., 2016); at this point, hypertrophy could occur if any of the remaining circulating SCs were still active.

Neutrophils

For initial time steps of the agent-based model, neutrophils were distributed randomly on ECM grid points. Once damage had occurred, neutrophils searched their neighborhoods for damaged ECM and moved towards these points by updating their location to points where collagen density had declined since the previous time step. Neutrophils peak within 24 h (Smith et al., 2008). When neutrophils encountered damaged ECM or fibers, they proliferated, released IL-6 (Xue et al., 2015), and marked objects as needing repair. Neutrophils then broke down nearby damaged objects.

Macrophages

At initial time points, macrophages were localized to ECM. Once damage had occurred in the model, macrophages searched a Moore neighborhood for the highest IL-6 concentration and moved towards the corresponding grid location. Once damage was located, the angle of movement was computed, and the macrophages moved in this angular direction towards the damage. This method simulates the behavior of chemical factors that attract macrophages (Serrano et al., 2008;Muñoz‐Cánoves et al., 2013). When a macrophage encounters a damaged fiber, the fiber is set to “needs repair” and is eligible for phagocytosis (Tidball, 1995; Oishi and Manabe., 2018). Phagocytosis allows for macrophage proliferation and increases in the levels of TNF-α present due to the proliferation of pro-inflammatory M1 macrophages in the model (Ostrowski et al., 1998; Xue et al., 2015; Wynn and Vannella, 2016). After clearing damage, M1 macrophages switch their phenotype to anti-inflammatory M2 macrophages which release TGF-β (Delaney et al., 2017; Dort et al., 2019).

Satellite Cells

In our simulations, SCs are seeded according to physiological values of approximately 0.10 SCs per fiber per 10 µm thick section (Verdijk et al., 2014; Snijders et al., 2015; Mackey et al., 2017). During initialization, SCs are randomly assigned to border fibrils to represent their location between the sarcolemma and basement membrane of muscle fibers. The activation of quiescent SCs requires the presence of hepatocyte growth factor (HGF) (O’Reilly et al., 2008). SC division can take place after 18 h for each agent (Zammit et al., 2002). The proliferation of active SCs then occurs to mimic the transit-amplifying cells present during repair (Hsu et al., 2014). This process took place in a timeframe when the gradient of IGF-1 was greater than zero (Zanou and Gailly, 2013). Both symmetric and asymmetric division occurred. Symmetric division resulted in two satellite cells or two myoblasts, while asymmetric division produced one quiescent and one active satellite cell (Kuang et al., 2008). The chance of division decreased 20% after each of the first three divisions for a single SC, before decreasing by 40% at the fourth division (Siegel et al., 2011). Active SCs searched Moore neighborhoods for fibers that need repair and placed a myoblast at these locations.

Fibroblasts

Fibroblast levels were seeded according to Mackey et al. (2017). Fibroblasts became activated myofibroblasts in the presence of a positive TGF-β gradient (Ismaeel et al., 2019). Myofibroblasts searched the area for empty cells that neighbored ECM components. Myofibroblasts competed with myoblasts to regenerate tissue by depositing collagen near damaged ECM edges. When there were no SCs present in the muscle, fibroblasts deposited collagen in any remaining spaces.

Secreted Factors

Levels of secreted factors IGF-1 (Martin et al., 2015), HGF (Leuning et al., 2018), IL-6, TNF-α (Kim et al., 2011; Martin et al., 2015; Xue et al., 2015) and TGF-β (Vignola et al., 1996) per hour were represented by the following equations:

where AM was the number of anti-inflammatory macrophages, N was the number of neutrophils, PM was the number of pro-inflammatory macrophages, Fb was the number of fibroblasts, and DE was the number of damaged ECM components.

Finite Element Modeling

The initial FE model was built from the same cross-sectional geometry that informed the initial agent-based model’s creation. Each subsequent FE model was built from the endpoint agent-based model geometry resulting from each ABM simulation. For the initial model, coordinate points from the segmented histological image were imported to Inventor^® (Autodesk., San Rafael, CA). The coordinates were then connected, and the muscle fiber bundle cross-section was extruded to represent 1 cm of muscle (Figure 1B). FE simulations were conducted in FEBio (Maas et al., 2012) and included muscle fibers and ECM. Muscle fibers were modeled with superposed active and passive stress, to simulate muscle activation in the fiber direction, where Cauchy stress is given by:

Muscle fiber stress, $T_a$, was modeled using the time-varying elastance model in FEBio:

Where tension of maximum isometric contraction $T_{max}=135.7kPa$, peak intracellular calcium concentration ${\left(C{a}_0\right)}_{max}=4.35\mu M$, $B=4.75\mu m^{-1}$, $l$ is the sarcomere length, and the length at which there is no active sarcomere tension is $l_0=1.58\mu m$. ECM and the passive component of the muscle fibres were modeled as nearly incompressible, hyperelastic materials, based on the following strain energy function for biological soft tissues (Weiss et al., 1996):

where $\mathbf{\Psi }$ is the strain energy functional, expressed uncoupled as the superposition of the ground substance Mooney-Rivlin response, ${\mathbf{F}}_1$, the response of the fiber family, ${\mathbf{F}}_2$, and the dilatational response where $K$ is the bulk modulus and $J$ is the Jacobian of the deformation tensor. The Mooney-Rivlin response is defined as:

where ${\tilde{\mathbf{I}}}_1$ and ${\tilde{\mathbf{I}}}_2$ are the first and second invariant of the deviatoric right Cauchy-Green deformation tensor, $c_1$ and $c_2$ are material parameters. The fiber response is

where $Ei$ denotes the exponential integral function, defined for real non-zero values below:

$\widetilde{\text{λ}}$ is deviatoric fiber stretch; and $c_3$, $c_4$, $c_{5 }$and $c_6$represent material constants of the constitutive relations. Parameters are provided in Table 1. Fiber stress response is defined as:

TABLE 1. Material properties for ECM and fibers in the FE models of the muscle fiber bundle. Properties were unchanged across all iterations.

Note that we used one set of constitutive parameters for both CP and TD FE models. Eccentric contraction was simulated by imposing a fiber stretch of 30% with one end of the muscle fiber bundle held fixed. Fibers underwent maximum voluntary active contraction during eccentric load. Activation parameters are given below (Table 2). At the end of FE simulations, element centers were calculated based on average coordinates for node points. Strain values were then matched with element centers before being mapped onto the 2D coordinate points of the agent-based model (Figure 3). Inhomogeneous strain values over the area of the fibers were parsed into the agent-based model.

TABLE 2. Active contraction properties.

FIGURE 3. Mapping FE element centers into ABM coordinate space. (A) FE element centers were calculated based on nodal points. (B) ABM pixels were imported with associated X, Y coordinates. (C) Element center coordinates and associated strain per element were then registered to ABM coordinates.

Coupled Model

Coupling of these models was achieved by simulating eccentric contraction and damage via strain using the FE models, followed by simulation of regeneration using ABM, where the endpoint geometry was then used for a new FE model. The strain values from the FE model were mapped to coordinate points of the ABM and the highest 10% of strain values were used to mark fibrils at relevant locations as damaged in the ABM following initialization (Figure 2). Two agent-based models were used for two separate simulation series: one agent-based model representing CP muscle and another representing TD muscle. Each agent-based model in our framework simulated 28 days of regeneration. The two seeded initial ABM environments represented TD muscle with an initial SC count of 10, or CP muscle with an initial SC count of 4. These simulation environments were run separately, and results from each agent-based model were used to create separate endpoint geometries: CP post-regeneration vs TD post-regeneration. Each endpoint geometry was exported to generate new FE meshes for the coupled simulation’s next iteration (Figure 4). FE constitutive properties were kept consistent across iterations. After each FE simulation, strain values were exported and parsed to the agent-based model, where they were registered with the correct coordinates (Figure 3). We chose a representative ABM geometry output near the mean, which was used to create the next FE model. This process was repeated three times to simulate 3 months of muscle damage and regeneration. Three months is the time it takes to investigate three full cycles of regeneration following injury, without the added complexity of increased simulation time and physiological processes such as growth. Since constitutive properties remained unchanged—across iterations and between TD and CP models—this approach probed the role that SCs have in muscle regeneration and the potential for chronic degeneration in the case of impaired SC density.

FIGURE 4. ABM and FE modeling coupling workflow. The initial pixel values and coordinates were used to build both starting ABM and FE modeling geometry. (A) FE simulations were run to model active eccentric contraction of a fiber bundle. (B) FE strain values were recorded with associated coordinate points. (C) Points of high strain were imported into the agent-based models as local damage. The agent-based models were run, and mean endpoint geometry (D) was used to generate a new FE geometry.

Results

Each agent-based model simulated 672 h (28 days) of muscle regeneration. Figure 5 illustrates the progression of regeneration following injury as simulated by our agent-based model. At t = 1, 10% of fibers were damaged, and neutrophil invasion began and lasted 24 h. This was followed by macrophage invasion over the first 3 days. Once macrophages cleared the necrotic tissue, SCs began the repair process. Fiber repair was completed at 216 h post-injury; hypertrophy occurred in cases where the environment was suitable, i.e., when there were still active SCs in the environment after the damaged fibers were replaced and boundaries re-formed. In these cases, border fibrils were re-formed by recalculating the final fibrils that were next to ECM, to accommodate hypertrophy. In the last stages of repair, fibroblasts acted to fill the remaining gaps with ECM.

FIGURE 5. ABM simulation of muscle regeneration over time. At the first time point, regions of the muscle fiber bundle show damage. This damage signals for the mobilization of macrophages, neutrophils, satellite cells and fibroblasts, which interact in the muscle regeneration process. Note that the dominant change from t = 216 to t = 672 is ECM remodeling.

FE material parameters ($c_1$, $c_3$, $c_4$, $c_5$, $K$ and density) were altered by ±10% to investigate changes in strain due to different input values for the material properties (Table 3). For both ECM and Fiber materials, a 10% increase in parameter $c_5$ decreased the maximum strain value by 7.07%; a 10% decrease in $c_5$ increased the maximum strain value by 7.29%. Percentage change in minimum strain was also altered by $c_5$ where a 10% increase in $c_5$ decreased the minimum strain by 2.16%, and a 10% decrease in $c_5$ increased the minimum strain by 2.17%. For parameters $c_1$, $c_3$, $c_4$, $K$ and density, ± 10% variations in input values resulted in less than 1% change in maximum and minimum strain values.

TABLE 3. Sensitivity analysis of FE material model parameters for ECM and Fibers. Model values were varied by ±10% (input values) and percentage changes to the maximum and minimum strain in each FE simulation are shown.

Damage levels with additional necrosis were varied from 0 to 20% at 5% increments to explore the range in which efficient muscle regeneration can be simulated (Figure 6A). At 5 and 10% damage levels, the clearance of damaged fibrils was completed within 90 h post injury and recovery of original fibril count was achieved. 5% damage led to endpoint fibril count of 10,020 ± 80 (mean ± SD) and 10% led to endpoint fibril count of 9,954 ± 206 (mean ± SD). 15% on average was able to clear the damage but unable to reach original fibril levels (endpoint mean ± SD 9785 ± 385). With extensive damage of 20%, there was insufficient clearance of damaged fibrils (Figure 6B, arrows), and repair that occurred around those damaged fibrils did not reach the original count of 9,864 fibers.

FIGURE 6. Fibril recovery post injury at varied damage levels of 0% fibrils to 20% with 5% steps. (A) Fibril recovery at 5 and 10% showed average end fibril count increase. At 15% clearance of damaged cells was effective and repair was close to recovery, however, at 20% damage, the clearance of damaged fibrils did not reach 20% (arrow) and therefore (B) inflammatory cells and damaged fibrils remain within the fibers.

We compared our simulated SC with literature values from Snijders et al. (2015) and Mackey et al. (2017), and inflammatory cell counts with Wosczyna and Rando (2018) (Figure 7). For SCs, initial and average peak values at 72 h are comparable; however, SC proliferation did not occur at the same initial rate compared to the Snijders et al. Simulated proliferation stalled for the first 24 h and then rose sharply while literature data showed a gradual rise in SCs. Simulated SC counts started to decline at 120 h, which was similar to the literature (Figure 7A). Average neutrophil counts in the agent-based models peaked at 20 h under typical ABM conditions, while macrophages peaked around 40 h. The decrease in macrophage count in the ABM simulations occurred in sync with the generalized time course for M1 macrophages, as evidenced from Wosczyna and Rando (2018) (Figures 7B,C). However, simulated macrophage count did not perfectly replicate the time course for M2 macrophages. This is partially due to the nature of the ABM, where the macrophage count combined both M1 and M2 macrophages as one cell population that underwent a switch at the appropriate time.

FIGURE 7. ABM simulated cell counts over time. (A) SC simulation results (mean ± SD) compared to Snijders et al. (Snijders et al., 2015) (B) Simulated inflammatory cell counts of neutrophil and macrophage (mean ± SD) over time (solid line) compared to the generalized time course of when neutrophil, M1 macrophage and M2 macrophage cell counts are above 20% from baseline (dashed line), adapted from (Wosczyna and Rando, 2018).

Fibril counts over the three iterations demonstrate the divergent outcomes between the TD and CP muscle milieu (Figure 8). In the first iteration ABM simulation for both CP and TD environments, the initial geometries were identical, and the initial fibril count was 9,864 for both models (Figure 8A). Damage was cleared at approx. t = 50 h and regeneration began at t = 70 h when SCs were activated. In the TD scenario, recovery over initial fibril count occurred at t = 200 h and a 0.79% increase in fibril count above original count (9,942 ± 108) was observed. In the CP seeded agent-based models, recovery peaked at 200 h but the model was not able to recover to original fibril count (9,296 ± 289). In the second iteration ABM simulation (Figure 8B), TD and CP environments were seeded with 10,000 and 9,515 fibrils, respectively, due to altered geometry that resulted from the first iteration. For both scenarios, a 10% damage threshold was applied again. Over the 28 days of simulated regeneration, the TD environment repaired its fibril count to just above original levels (0.02%; fibril count = 10,002 ± 130). Over the same 28 days, the CP environment did not repair to original levels of fibrils and reached a final fibril count of 9,081 ± 331. The initial fibril counts for the third iteration (Figure 8C) of the ABM simulations were 10,060 (TD) and 9,272 (CP). For both cases, repair began at t = 78 h. TD fibril count reached a fibril count of 10,080 ± 260, which is above the initial fibril count for the third iteration. CP recovery was again unable to repair to original counts, completing its regeneration cycle at a fibril count of 8,763 ± 301.

FIGURE 8. ABM cell counts over time (mean ± SD) seeded with CP (SC = 4) or TD (SC = 10) initial conditions. (A) Iteration one cell counts over the first 672 h, based on initial geometry. TD fibril count exceeded initial count whereas CP fibril recovery was impaired. (B) ABM cell counts (mean ± SD) using TD and CP iteration two geometry and strain values seeded with CP (SC = 4) or TD (SC = 10) conditions. TD fibril recovery continued to exceed initial counts and CP fibril recovery was further impaired. (C) Iteration three of the ABM cell counts. TD fibril count peaked above the original value during repair however stabilised to just below initial values by the end of the simulation. CP fibril recovery decreased to 8,763 ± 315 (mean ± SD).

The coupled agent-based model and FE model simulated cyclic damage and regeneration in TD and CP muscle, which manifested as progressive degeneration in the CP model. When the FE model geometry simulated eccentric lengthening of muscle, the regions of highest strain occurred on outer fibers, particularly those located on the corners of the geometry. These regions were then assigned high strain and damage in the agent-based model. For the TD coupled model, i.e., the model in which the simulated muscle milieu represented TD muscle, the agent-based model repaired all of the damage and grew larger than the initial geometry with an increase in mean muscle fraction from 84.1 to 84.7% (p<<0.001), demonstrating emergent hypertrophy in these simulations. For the CP coupled model, i.e. the model in which the simulated muscle milieu represented CP muscle, the agent-based models were unable to repair the damage induced by the FE model (Figure 9) and the mean muscle fraction declined from 84.1 to 79.2% (p<<0.001). In the second TD iteration, the FE model displayed more widely distributed high strain values on outer fibers compared to the first iteration. In this second iteration of the TD model, the agent-based models again showed an increase in size of fibers affected by damage (muscle fraction = 85.2%, p<<0.001), again simulating emergent hypertrophy. In contrast, in the second iteration of the CP model, the highest strain in the FE model was concentrated on the smallest outer fibers, and in the ABM simulations, this damage was not entirely repaired, failing to restore the geometry to its original geometry and further decreasing muscle fraction to 77.3%. In the third iteration, fiber geometry showed a marginal increase in size in the TD model. The CP fibers continued to decrease in size in the third iteration, where muscle fraction fell to 79.6% and three fibers in particular demonstrated considerable atrophy (Figure 9 arrows). Average fiber cross sectional area increased to 3,251 µm² for TD scenario, and decreased to 2,533 µm² for CP, from original area of 3,182 µm², over 3 months.

FIGURE 9. In the coupled ABM-FE modeling mechanobiological simulations, the initial geometry leads to two unique endpoint geometries, according to whether the agent-based models had a cellular milieu based on CP or TD muscle. Each of those endpoint geometries then leads to a new FE model geometry and simulation, where high strains differ based on the geometry from the previous step. Ultimately, divergent geometries emerge, reflecting the different CP vs TD muscle outcomes, where TD muscle regenerates fully each cycle and CP muscle fibers cyclically degenerate.

Sensitivity analysis for SC counts of 4 (CP), 5, 7, 10 (TD) and 13 were performed and evaluated with respect to fibril recovery (Figure 10). Mean fibril recovery increased with seeded SC count. SC counts of 4, five and seven had greater variance compared to SC counts of 10 and 13. Higher seeded SC count was inversely related to the number of hours required for recovery of initial fibril count. A 50% reduction of SCs reduced average fibril recovery count by only 110 fibrils while a 30% reduction in seeded SC count led to mere marginal decreases in endpoint fibril recovery. An SC count of 4, however, resulted in a more dramatic decrease of 523 fibrils compared to an SC count of 5.

FIGURE 10. Sensitivity analysis for SC count and the effect on fibril count (mean ± SD) over simulation time course. SC count was set to 4 (CP), 5, 7, 10 (TD), and 13. Mean fibril count increased with increase in seeded SC count.

Endpoint ECM count was measured after the third month of regeneration for CP and TD environments (Figure 11A). The CP environment had an endpoint ECM count of 2,240 ± 176 and the TD environment had a lower endpoint count of 1,932 ± 139. Tissue composition was from the initial simulation was calculated and compared to both the third month averages for CP and TD simulations (Figure 11B). The third month endpoint CP ECM made up 20.4% of the tissue, increased from 15.8%. In comparison, the third month endpoint TD ECM increased marginally to 16.1% from initial. While the TD ECM fraction and muscle fractions were stable over 3 months, the CP muscle fraction decreased from 84.2% to 79.6%.

FIGURE 11. (A) ECM recovery over the third month 28 days post injury simulation. CP muscle environment (SC = 4) had higher end point ECM count of 2,240 ± 176 and the TD environment had a lower endpoint count of 1932 ± 139. (B) Tissue composition changed from initial simulations to the end of the third month. ECM percentage increased by 4.6% in the CP scenario and marginally increased by 0.3% in the TD scenario.

Discussion

In this work, we coupled agent-based modeling with finite element modeling of a muscle fiber bundle to explore the interaction of muscle damage and regeneration in the context of altered satellite cells and supporting cells in the muscle fiber environment. We were particularly interested in whether this framework could demonstrate a progressive degeneration of the muscle fiber bundle consistent with the progression of cerebral palsy.

Our coupled model demonstrates the canonical process of muscle regeneration after several bouts of damage from eccentric contraction (Charge and Rudnicki, 2004; Toumi et al., 2006; Järvinen et al., 2008; Wang and Rudnicki, 2012; Novak and Koh, 2013). Over the course of 672 h, our simulations showed damage, inflammation, clearance of damaged tissue, repair of muscle tissue regions, and remodeling of the extracellular matrix. The regeneration model includes quiescence, activation, and proliferation of satellite cells; M1 and M2 macrophages and neutrophils; fibroblasts; muscle fibers and extracellular matrix agents; and secreted factors of TNF-α, TNF-β, IGF-1, and HGF. Each of these factors and agents interacts with one another in a bottom-up and nonlinear approach. The timeframe for repair in the typical scenario was 14–21 days, in line with active muscle regeneration studies where repair peaks at 2 weeks and subsequently declines (Ambrosio et al., 2009). Here, we imposed 10% damage at the initiation of each simulated iteration. It is possible that more extensive damage would have required a longer simulation time frame to repair fully; however, consistent in the literature are specific timeframes for different stages of repair. This suggests that the cellular response may scale in magnitude with the level of damage, rather than with the response timeframe. With this said, there are certainly practical limits on the extent of damage that can be repaired in skeletal muscle that is shown by incomplete clearance of damaged tissue in the 20% damage scenario. Prolonged presence of damaged fibers and resultant cellular debris due to insufficient inflammatory cell response has been shown to delay muscle regeneration (Summan et al., 2006). Our choice of 10% damage is thought to represent a reasonable level of mechanical damage that does not exceed skeletal muscle’s typical ability to mount the repair process.

Our models were consistent with literature reports of peak satellite cell concentrations and timing of inflammatory cell concentrations (Snijders et al., 2015; Wosczyna and Rando, 2018), and demonstrate other well-characterized phenomena observed in skeletal muscle regeneration such as ghost fibers (Webster et al., 2016) and fibrosis (Booth et al., 2001; Peterson et al., 2012; Kinney et al., 2017). The rules for the present model focus on post-injury cell counts during regeneration in healthy individuals after exercise-induced mechanical damage (McKay et al., 2010; Mackey and Kjaer, 2016; Kim and Lee, 2017; Nederveen et al., 2017). While this is a suitable model for exploring typical muscle injury and regeneration, this remains a largely unexplored area in the context of CP. Additionally, directional cues for non-cellular systemic agents such as cytokines have not been well-described in the literature. Lacking this information, the present models relied on cytokines that were concentrated at damage sites. Future models might consider using geometrical features to inform directional cues for the diffusion of signaling factors. For example, the inclusion of capillary location within the geometry could enhance the spatial localization of satellite cells and systemic cytokines within the model. However, at this stage, the functions of systemic and local signaling factors on muscle repair are yet to be fully explained (Chargé and Rudnicki, 2009); future experimental exploration and inclusion in computational models will be an exciting area in this field.

Sensitivity analysis demonstrated relative insensitivity of fibril count to satellite cell concentration between SC levels of 5 and 13; however, the level of five was such that, on average, simulations did not fully repair the damage induced in them and fibril count was lower at the end of the simulation. SC counts below five had drastically impaired ability to repair damage. The behavior of satellite cells in our model suggests a threshold of satellite cell concentration, above which is sufficient to repair 10% damage, and below which the damage repair falls off sharply. We performed additional sensitivity analyses on the constitutive parameters in our finite element model. While strain values were relatively insensitive to perturbations in most constitutive parameters, predictably we found a larger dependence of strain on c5, the along-fiber modulus. Overall, since the role of the finite element model was to identify regions of high strain within the geometry, the exact values of constitutive parameters are not thought to be of central importance to this goal. In this coupled modeling framework, we used strain as the mechanical parameter associated with damage in the muscle fiber bundle. Several prior studies reinforce this association (Lieber and Friden, 1993; Best et al., 1995; Lieber and Fridén, 1999). Mechanobiology is an active field, however, and it is unclear how other parameters such as strain rate, repeated loading, shear, and stress may contribute to the signaling and processes of muscle regeneration. Future modeling and experimental work may explore this area in finer detail.

Our framework involving FE modeling and ABM demonstrated a progressive degeneration of muscle fibers in the simulations resembling cerebral palsy muscle. The cerebral palsy muscle here was denoted by a reduced satellite cell concentration. Here, for simplicity and direct comparison, we did not alter material properties between the CP and TD FE models, and we did not change material properties over the course of simulations. It is arguable that CP muscle would have different constitutive properties compared to TD, such as stiffer muscle, or that changes over time to the muscle could well be captured with changes to material properties. The mode of injury used in this model was also limited to eccentric lengthening. Concentric exercise can result in hypertrophy; however, eccentric exercise is more efficient at eliciting hypertrophy and muscle growth (Schoenfeld et al., 2017). In any case, given the same form and extent of injury, whether by concentric or eccentric loading, the repair mechanism would likely remain the same.

In this model, we chose eccentric contraction as a standard model that causes fiber damage (and initiates the muscle repair process) and is physiologically relevant. From that stimulus, we simulate and observe the regeneration process. Damage frequency from eccentric strain may be caused more than once per month, however this time course was used to explore the entire four step process of muscle regeneration that is commonly referred to in the literature. Additionally, it is known that strain is localised during eccentric contraction and subsequent injury occurs in the areas of highest strain (above a threshold) (Lieber and Friden, 1993; Best et al., 1995; Garrett, 1996). In a test case, random damage was seeded across the same agent-based model (Supplementary Figure S1). This model was unable to clear cellular debris at 10% damage level and prevented ECM remodeling. This suggests that muscle repair in response to eccentric contraction vs inflammatory myopathies require different responses, the latter of which is beyond the scope of this study.

The purpose of the current model was to demonstrate if the process of muscle regeneration could lead to degeneration with reduced SC numbers. Under this proposition alone, our model illustrated this phenomenon, and the endpoint geometries observed here are similar to those seen in the literature for cerebral palsy muscle cross-sections (Marbini et al., 2002; Foran et al., 2005; De Bruin et al., 2014), namely reduced muscle fiber cross-sections and increased area fraction of ECM. To make a more robust geometry selection for the subsequent iterations, an average pixel geometry could be generated in the future versions of this model. Over 3 months, the muscle fraction declined in cerebral palsy simulations by 4.6%. In these simulations, the change in muscle geometry due to decrease in fiber size of damaged fibers represents atrophy at the fiber level. Under the framework that we have presented here, this degenerated muscle architecture characteristic of CP emerges from the same initial muscle geometry as in the TD model. The only difference was the muscle milieu, i.e. number of satellite cells, and thus the regeneration process of the muscle after injury.

Conclusion

Coupled modeling is a powerful tool in its ability to connect tissue and organ level behaviors to simple cellular interactions. This model demonstrated growth over time in a TD muscle environment that experienced strain similar to that which would occur from active eccentric lengthening. In simulations of CP muscle environment, the same strains led to gradual degradation of size and shape of muscle fibers over time. Overall, this work suggests a plausible connection consistent with the physiological mechanisms that are observed in the clinical manifestation of cerebral palsy.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

SK developed and analyzed all models presented in this work, is the author of the code associated with model development and analysis, and was the principal author of this manuscript; JF collaborated on modeling approaches, overall conceptualization of this work, and contributed insight and critical review; GH is the corresponding and supervising author, developed initial concept, oversaw modeling, and contributed insight and critical review.

Funding

Funding for this work was generously provided by the Robertson Foundation Aotearoa Fellowship Award Number 3715249.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2021.689714/full#supplementary-material

References