The FRiND Model: A Mathematical Model for Representing Macrophage Plasticity in Muscular Dystrophy Pathogenesis

Muscular dystrophy describes generalized progressive muscular weakness due to the wasting of muscle fibers. The progression of the disease is affected by known immunological and mechanical factors, and possibly other unknown mechanisms. This article introduces a new mathematical model, the FRiND model, to further elucidate these known immunological actions. We will perform stability and sensitivity analyses on this model. The models time course results will be verified by biological studies in the literature. This model could be the foundation for further understanding of immunological muscle repair.


Introduction
Muscular dystrophy (MD) is a group of genetic disorders whose archetypal pathology is progressive weakening of skeletal muscle tissue. Severity of degeneration and muscles affected vary depending on the type of MD (and person) with some forms leading to early death while other forms remain unnoticed until adulthood. All forms of MD combined affect about 37 per 100, 000 individuals in Northern England (Norwood

Concept of the FRiND Model
To understand the interactions needed for muscle repair, we could create a model that incorporates three groups of cells: macrophages (in number of cells per mm 3 ), muscle tissue cells (in percentage of tissue), and myocytes (in number of cells per mm 3 ). The interactions of these cells in a mathematical model should ideally stay in equilibrium until any damage to muscle occurs (Figs. 1,2,3). M i is activating a macrophage phenotype.

Fig. 1
Macrophage interactions This figure describes the interaction of different types of macrophages during muscle repair. Inflammatory signals from resident macrophages are released after tissue is damaged prompting invasion by pro-inflammatory macrophages (M 1 ). Phagocytosis of damaged tissue by pro-inflammatory macrophages promotes inhibiting macrophage phenotype (M 2c ) which starts deactivating (M 3 ) pro-inflammatory macrophages. Resolution macrophages (M 2a ) appear after further phagocytosis of damaged tissue by inhibiting macrophages. M 2a macrophages slowly apoptose during muscle repair and following resolution of repair will return to starting population of resident macrophages. Conflicting signals by pro-and anti-inflammatory macrophages force some macrophages into a phenotypeless state (M 3 ) For purposes of this model, macrophage populations have been split into four phenotypes derived from in vitro studies: M 1 , M 2c , M 2a , and M 3 . Although in vivo studies have shown few macrophages display the exact qualities of these phenotypes, they can be shown to have certain genes-corresponding to those in vitro phenotypesupregulated (Novak et al. 2014).
Since monocytes infiltrate into damage tissue as pro-inflammatory macrophages (M 1 ,) the M 1 population increases proportionally to the damage dealt by injury. M 1 macrophages release cytokines T N F α and IL-1, which further allow M 1 infiltration and recruitment of phenotypeless macrophages (M 3 ); this is modeled using k 1 D 2 . These cytokines also promote pro-inflammatory genes in resident macrophages forcing them into an M 3 phenotypeless state. TGF-β is released when macrophages phagocytose damaged tissue and apoptotic immune cells; this promotes a shift in macrophages to an anti-inflammatory phenotype (Arnold et al. 2007).
Anti-inflammatory macrophages (M 2 ) fall into two phenotypes: the deactivating macrophages (M 2c ) and resolving macrophages (M 2a ). Increased exposure to TGF-β and certain glucocorticoids promote M 2c phenotype genes in M 1 macrophages. IL-10 which is released by M 2c macrophages inhibits the effects of IL-1 and suppresses M 1 macrophage's inflammatory abilities effectively becoming phenotypeless. Further exposure by M 2c to TGF-β promotes a shift to the M 2a phenotype. M 2a macrophages release cytokines IL-4 and IL-13 which accelerate M 2c phenotype change and promote myoblasts differentiation (Zhang et al. 2016).
Phenotypeless or "switch" macrophages (M 3 ) are used in this model to represent any macrophage which has a phenotype not sufficiently polarized or is suppressed by  another macrophage. The M 3 macrophages will act only as a feedback loop for M 1 and M 2a (Malyshev and Malyshev 2015;Kalish et al. 2017). Following damage to normal muscle tissue (N ), damage tissue (D) increases. Macrophages phagocytose this tissue and leave the area ready for regeneration (R). Transforming Growth Factor Beta, TGF-β, which is released by M 2 macrophages promotes the formation of fibrous tissue (F) (Arnold et al. 2007;Mounier et al. 2013). This tissue is the foundation upon which myotubes bind to form new healthy muscle tissue (N ) (Ogawa et al. 2015). In the model, the F, R, N , and D are given as percents of overall tissue as this is assumed to be conserved.
A subpopulation of satellite cells called myocytes form the core of new muscles. Pro-inflammatory signals from M 1 macrophages activate MyoD in resident satellite cells becoming myoblasts allowing for proliferation. The model simplifies this by combining the resident satellite cells and myoblasts into M b . When exposed to antiinflammatory signals by M 2a , myoblasts lose the expression of Pax and differentiate. These differentiated myoblasts later form into myotubes. The model combines the differentiated myoblasts and myotubes as a single variable, M d . M 1 macrophages are known to signal M d myoblasts to prematurely apoptose (Ogawa et al. 2015).
Damage to muscle tissue is controlled by a time-dependent function f (t) which can be designed to simulate different types of muscle injury. For acute injury, f (t) could be set as a Gaussian normal or a Gaussian function as described by Jarrah et al. (2014). For smaller but periodic injury, f (t) could be set with a random amplitude and period sin 2 function. This last function will be discussed in Sect. 3.3.
We will call this model the FRiND (Fibrous, Regenerating, inflammation, N ormal, Damaged) model.

Equations
The FRiND model uses the law of generalized mass action (GMA) to convert the diagrams in Sect. 2.1 into ordinary differential equations (ODEs). GMA states that the rate of reaction/interaction is directly and/or indirectly proportional to the product of the numbers (or masses) of objections. For each interaction, a kinetic parameter is incorporated and is estimated using data. In the FRiND model, activation and deactivation will be modeled using direct proportions; inhibition of an interaction will be modeled using an indirect proportion. An extra term −g(t)M 1 has been added for predictions of Diphtheria toxin in Sect. 3.2; we assume g(t) = 0 in all other sections.

Estimating the Model Parameters
The FRiND Model has 20 parameters plus any parameters needed for the initial damage function. Data for estimating the parameters come from Mounier et al. (2013) and Ogawa et al. (2015) (see Table 1). The choice of these literature sources was influenced by the similarity of the data collected with respect to the model proposed here; other sources include the adaptive immune system. Both experiments used cardiotoxin (CTX) as a source of damage to mouse muscle and similar control mice. The damage by CTX is modeled using a Gaussian function: A Gaussian function was chosen to model CTX to allow for a symmetric growth and decay with a predictable peak within one day that causes damaged tissue to reach over 80% (Mounier et al. 2013;Hardy et al. 2016). However, nearly any bell-shaped or log-bell function in time works when applied to the FRiND model. The data were collected from a small number of figures pertaining to the control mice in each paper. M 3 macrophage numbers were developed by subtracting the total macrophage populations by the number of M 1 and M 2 macrophages. To make estimating parameters easier, interpolation was used to calculate myocyte data on days when tissue data were present. The data for days 14 and 24 were the same as initial conditions since both papers report resolution of repair within 14 days. A standard genetic algorithm (Bäck et al. 1997) and pattern search (Hooke and Jeeves 1961) were used to train parameters to the data on both COPASI (Hoops et al. 2006) and MATLAB (The MathWorks, Inc., Natick, MA, USA).
The parameters d 3 , d 4 , and d 5 are misleading when first observed. These values represent M 1 , M 2c , and M 2a macrophages, respectively, ability to perform phagocytosis not the cytotoxicity. In Mounier et al. (2013), they demonstrated that in AMPKα1 − /−, mice phagocytosis of damaged tissue was hampered, and M 1 macrophages did not transition to M 2 macrophages. This leaves the possibility that M 2 macrophages had a greater capacity to phagocytose. This has also been speculated in phagocytosis of tumor cells by M 2c (Herter et al. 2014) and phagocytosis by M 2a in neurological repair (Ghosh et al. 2016). Multiple tries were made to estimate the parameters with d 3 > d 4 , and all failed to find a satisfactorily low objective value.

Fixed Points
A steady state is a multi-dimensional point where the system will remain constant over time. This occurs when the pointx (also called a fixed point) gives dx dt = 0 for all x. For a large dynamical system, finding a fixed point analytically can be difficult (Strogatz 2018). In the FRiND model, several simplifying assumptions can be made.
(1) Since the change in substrate amount is given in the model, parameters may be assumed to be always positive and nonzero. Furthermore, (2) steady state and initial conditions of the substrate must be nonnegative to make biological sense. (3) We will, though, insist that M b (0) > 0 (the initial condition of M b myocytes) as this will allow satellite cells to proliferate. Finally, (4) we want to see how the system responds after contractual damage has dissipated or any diphtheria toxins have left; this would require lim t→∞ f (t) = 0 and lim t→∞ g(t) = 0. Although the initial conditions for substrates depend on the situation being modeled, the muscle cells have a built-in conservation law dF dt Theorem 1 Assume the FRiND model follows the above assumptions and the parameters given in Table 2, then the model will have these two steady states.
Proof Using the definition of fixed points, we have: Consider Eq. (15) eitherD = 0 or the macrophage populations would be zero. If M 1 = 0,M 2a = 0, andM 2c = 0, then by Eq. (16) we would haveD = 0. Summing Eqs. (16)-(19) and usingD = 0, we have: Hence, we have a contradiction. This means thatM 1 ,M 2a , andM 2c cannot all be zero. Notice that Eq. (22) 12) and (14) are linearly dependent and give us a free variable,F. This give us: Assume now thatF = 0, then by Eq. (16) we have that eitherM . Thus,M d = 0 by Eq.(21). Notice, however, this gives you a fixed point of which is a special case of the earlier fixed point.
IfM 3 = k 3 k 4 M 2c (0), then substituting into Eq. (18) we have: However, by Eq. (20), this leads to: Using Eq. (21), we have: . Therefore, we have that our fixed point as In the fixed pointx 1 , we have a free variableF. This indicates that the associated steady state depends on when F(t) (the transient function of fibrous tissue) becomes constant. Since these values when near steady state-notice that regenerating tissue should already be close to zero-depend only on M d (t) (the transient function of differentiated myocytes), the FRiND model shows that the replacement of normal tissue by fibrous tissue depends on locally available differentiated myocytes. This observation is in agreement with Ogawa et al. (2015) and the model created by Virgilio et al. (2018). Normally, this is not an issue since myoblasts proliferate enough during the pro-inflammatory stage to give enough differentiated myocytes to repair muscle. However, this process can be disrupted enough for the free variable to become an issue as we will see in Result Section. Notice thatx 2 can violate assumption (2). This would allow a model to approach a steady state that is unrealistic. Stability analysis will show that the FRiND model will generally repel from this outcome. Theorem 2 will give a handy criteria forx 2 existing.
Therefore,M b exists and is nonnegative. Furthermore, we have that Hence,M d exists and is nonnegative.
Therefore, all components of the second fixed point exist and are nonnegative. The FRiND model will have a second steady state.

Stability Analysis
A steady state can be called stable, unstable, or saddle depending on the system's reaction to small perturbations to the fixed point (other types of reactions can occur but will not be needed in this discussion). The system will return to a steady state after perturbations when a fixed point is stable, whereas the system will repel if the fixed point is unstable. With a saddle fixed point, the system will repel from the point unless on specific vectors (Strogatz 2018).
The sufficient conditions for stability are given by eigenvalues of the Jacobian evaluated at the fixed points. A stable steady state will have only real, nonpositive eigenvalues. An unstable state will have only real, nonnegative eigenvalues, and a saddle steady state will have a mix of positive and negative, real eigenvalues (Strogatz 2018).
The eigenvalues of the Jacobian evaluated atx 1 are According to the assumptions listed earlier, the parameters and initial values are positive. Thus, the eigenvalues are real and nonpositive. Therefore, the steady state atx 1 is stable.
The eigenvalues of the Jacobian evaluated atx 2 are Again by the earlier assumptions, M 2c (0), a 1 , k 3 , k 6 , k 4 , and k 7 are positive. Hence, all the eigenvalues are real, and the eigenvalue M 2c (0)a 1 k 3 k 6 k 4 k 7 is always positive and −k 9 is always negative. Therefore, the steady state atx 2 is a saddle node.
The above stability analysis shows that the FRiND model will converge tox 1 unless on specific vectors (which will converge tox 2 ). The flexibility given by the free variable inx 1 will allow for the comparison of the normal steady state for the immune system to steady states that include replacement of normal muscle fiber with fibrous connective tissue.

Identifiability and Sensitivity Analysis
The question could arise at this point whether the FRiND model's parameters can be estimated by any appropriate dataset. In mathematics, this concept is called identifiability. Models which are unidentifiable can be reduced by eliminating parameters; this helps to estimate parameters faster and creates a uniqueness to parameter estimation. This concept is related to sensitivity analysis. A system which is highly sensitive will be at least locally identifiable (Eisenberg and Hayashi 2014;Fisher 1959).
To calculate the sensitivity matrix numerically, we estimated ∂ x k ∂ p j (t i ) for each i, j, k with the central-difference formula. We also tried more accurate estimations of the partial derivatives before writing this article but found the results to be similar. We calculated the sensitivity matrix and Fisher Information Matrix (FIM) as described in Eisenberg and Hayashi (2014 This gives us the determinant of the FIM matrix as 1.0422 · 10 226 . As the determinant is large, the parameters of the FRiND model are sensitive and at least locally identifiable (Eisenberg and Hayashi 2014;Fisher 1959).

Results
One of the goals of the FRiND model is to use a damage function, f (t), to study the effects of damage to muscle tissue from chronic immune activation; specifically, the FRiND model allows exploration into how the acute immune response becomes a chronic immune response which causes muscle to be replaced by fibrous tissue.
To accomplish this goal, we must show that that model is accurate on the given data used to estimate the model parameters (Sect. 2.3) which were obtained from an acute damage event. After the accuracy of the model for given data has been shown, we will use the FRiND model to predict other experiments and results for which the model has not been trained including simulations of chronic muscle damage and effects of diphtheria toxin on muscle repair. All time-course graphs were created using Adams-Bashforth-Moulton Predictor-Corrector method (Cheney and Kincaid 2007) on MATLAB and verified using COPASI (which uses LSODA) (Hoops et al. 2006). Table 1 As mentioned above, parameters were estimated using data from the literature (Mounier et al. 2013;Ogawa et al. 2015) (Table 1). The parameter values are reported in Table 2. As discussed in Sect. 2.3, we chose the Gaussian function (Eq. (11)) to simulate damage from Cardiotoxin (CTX) (Dell'Acqua and Castiglione 2009;Jarrah et al. 2014).

The FRiND Model is Consistent with the Data in
Muscle tissue predictions followed closely to that given by Mounier et al. (2013). From Fig. 4a, normal tissue experiences a steep drop as muscle is exposed to CTX; normal tissue is completely compromised within 24 h. Damage peaks within the first 24 h followed by complete phagocytosis over the next 3 days. The model gives regenerating tissue at about 40% at day 2 which is slightly less than phagocyted myofibers given in the literature (Mounier et al. 2013). However, the model's values for regenerating tissue at days 4 and 7 are closer. Fibrous tissue occurs earlier and peaks slightly lower than that reported by the literature (Mounier et al. 2013). Recovery of new normal tissue does occur within the reported period (Fig. 5). The model also closely follows the trends in tissue data reported by Arnold et al. (2007) which used a notexin (Fig. 6).
Macrophage predictions followed the general characteristic of the data given in Mounier et al. (2013). From Fig. 4b M 1 , macrophages peak between day 1 and day 2 as reported by Mounier et al. (2013); Arnold et al. (2007). However, the number of M 1 macrophages is less than expected. M 2 macrophages are noted to peak between days 2 and 4 as expected (Mounier et al. 2013;Arnold et al. 2007); the numbers do, however, understate day 2 and overshoot day 4 numbers. M 3 macrophage predictions closely follow data from Mounier et al. (2013).
Myocyte predictions closely follow the expected cell numbers before day 7 (Fig. 4c). After day 7, myoblast numbers fall short of expected numbers and differentiated myoblast numbers overshoot reported cell numbers (Ogawa et al. 2015).

The FRiND Model is Consistent with Inflammatory Macrophage Inhibition by
Diphtheria Toxin Arnold et al. (2007) performed several experiments studying the actions of macrophages on muscle regeneration. The group noticed when M 1 macrophages are depleted within the first day after damage by the exposure to diphtheria toxin (DT) that necrotic (damaged) tissue removal is delayed significantly. Furthermore, the regeneration process resumes after M 1 macrophages were allowed to return and the steady state lacked significant replacement of normal tissue by fibrous tissue. To simulate this, the FRiND model can be modified by subtracting g(t)M 1 from the M 1 macrophage equation. Thus, where g(t) = 50 for 0.3 < t < 2 0 for otherwise . Figure 6 shows the results of this adapted FRiND model. M 1 macrophages are depleted shortly after diphtheria toxin is administered and rebounds around day 2 when k(t) = 0. While M 1 macrophages are depleted, only a small amount of necrotic tissue is removed mainly by M 2 macrophages already in the tissue. The repair process resumes after DT subsides and M 1 macrophage numbers rebounds with full recovery by 14 days. The steady state of the system also shows that no significant fibrosis remains by 14 days. This matches the results given by Arnold et al. (2007).

The FRiND Model is Consistent with Chronic Immune Activation and Literature Results in mdx Mice
The archetypal pathology of MD is the replacement of healthy muscle tissue with fibrous tissue caused by severe acute damage to muscle from normal daily activity (Desguerre et al. 2009). Most studies of MD and muscle repair revolve around performing a single acute muscle damage. However, MD patients/mice rarely suffer from just a single damage event; instead normal daily activities cause both mild and severe acute events throughout the body. We would like to use the FRiND model to predict this condition by only changing the input damage function. This could show that repeated severe acute damage events cause the chronic immune response and replacement of muscle tissue, which precipitates the muscle weakness in MD.
A choice for this daily damage could be f (t) = a sin (bt + c) for some choice of parameters a, b, c. However, to avoid negative damage we should use f (t) = a sin 2 (bt + c). Since we can dictate whether t = 0 is at midnight or midday, without loss of generality, we can choose b = π and c = 0 which would give daily periodic damage peaking at midday or midnight. To simulate normal daily muscle usage, this model will allow a to be a random number between 0 and 1 which will be selected each day. For Theorem 1 to apply in the FRiND model, damage added by f must either converge to zero naturally or by setting the function to zero at some time point.
The following results are again given using Adams-Bashforth-Moulton Predictor-Corrector method on MATLAB and f (t) = a sin 2 (π t) . Figure 7 shows that over a period of several weeks, normal tissue is repaired with some replacement by fibrous tissue. M 1 macrophages are nullified by M 2c macrophages and prevented from encouraging proliferation of more myoblasts. This causes a lack of differentiated myoblasts to bind with fibrous tissue into normal tissue; a buildup of fibrous tissue follows. Normal tissue is partially replaced by fibrous tissue around 16 weeks. A partial replacement by fibrous tissue matches studies of gastrocnemius muscle in mdx/utrn −/− mice (Lu et al. 2014) which found around 20% of normal tissue was replaced by collagen-positive tissue after 8 weeks.
The phase portrait of the FRiND model ( Fig. 8) when f (t) = a sin 2 (π t) emphasizes the cyclic nature of the process. Fibrous and normal tissue interchange constantly creating jagged ellipses throughout the first 11 weeks. Following the depletion of differentiated myoblasts, the phase portrait shows that the cyclic nature ends and spirals toward the F(t) = N (t) − 100 line which is a push toward a steady state with elevated fibrous tissue and loss of normal tissue (both damage and regenerating tissue will approach zero by the 16th week. A similar function to consider is f (t) = a sin 2 (π t) t+1 ; this function allows for weakening of a patient's strength over time as muscle is replaced by fibrous tissue. Using this damage function, we can create results (Figs. 9 and 10) with similar properties as the earlier non-fading damage function and fibrous tissue reaching steady state around 7% which closely aligns with mdx mice studies which found 4%/7% (male/female) (Salimena et al. 2000) and 7% (Sun et al. 2015). few mathematical models have been created to quantify the experiments about MD (Houston et al. 2018).
Although most forms of MD are genetic disorders, chronic inflammation caused by repeated severe acute immune responses has been implicated in several forms of MD (Wehling et al. 2001;Selva-O'Callaghan et al. 2006;Spencer et al. 2001;Desguerre et al. 2009;Tidball and Villalta 2010). The two previous models (Dell 'Acqua and Castiglione 2009;Jarrah et al. 2014) attempted to simulate research of inflammation during the lifetime of mdx mice (Wehling et al. 2001;Spencer et al. 2001). Both models used the law of mass action with three immune cells: CD8+ and CD4+ T cells and pro-inflammatory macrophages. They used a single acute damage event to start a perturbation that estimated chronic inflammation. The models were also bistable; either the normal tissue completely returned or a little damage remained permanently.
However, the previous models miss several key aspects of tissue repair. Treating macrophages as a homogeneous group reduces the predictive power of those models. Macrophages display several phenotypes that perform different tasks during muscle repair; some of which play an important role in laying fibrous connective tissue and migrating nascent myotubes to form new muscle; previous models handled this by using linear mass action dynamics for tissue repair after damaged tissue was eaten by pro-inflammatory macrophages. We introduced the FRiND model to rectify these issues.
The FRiND model offers many avenues for future usage by integrating macrophage plasticity and nonlinear mass action dynamics. The ability to adapt to experiments with different criteria for subject injury and to study different pharmaceutical outcomes by changing parameters in the system could be useful in exploring biological outcomes in muscle regeneration.
We have shown that changing f (t) in the FRiND model can be used to explore outcomes that depend on the injury caused to muscle tissue. In addition, this injury can be different than the injury used to give the estimating parameters. The parameters used throughout the results section were estimated using control mice after CTX damage. The results studied, however, were free to change the f (t) to explore other forms of muscle damage like normal degradation due to muscular dystrophy. This could hint that the weakening of skeletal muscle in several types of MD is the result of the immune system's overreaction due to compromised muscle cell structure.
The FRiND model also predicts experimental and pharmaceutical outcomes by changing or adding parameters into the model. As shown in Sect. 3.2, an added parameter in the model exhibits the same general trends reported by Arnold et al. (2007). This could allow future outcomes to be tested and studied before performing an experiment.
The FRiND model does include some simplifications. The model does not include several immune and non-immune cells which have been shown to be significant in tissue repair including neutrophils (Arnold et al. 2007), T cells (Tidball and Villalta 2010;Burzyn et al. 2013), fibro/adipogenic progenitors (Joe et al. 2010), and tissue cells under inflammation. With more data, many of these cells could be incorporated into the FRiND model. Inclusion of these cells and many others could allow further elucidation into the repair process.
The FRiND model also simplifies spatial dynamics by following the law of generalized mass action which assumes interactions are isotropic. Future research could expand this system into an anisotropic model using partial differential equations (PDEs). This could also allow a three-dimensional study of muscle regeneration. Another direction could include stochastic differential equations to differentiate early repair when cells are first infiltrating from later stages when cells are nearly isotropic.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.