Juan B.Gutierrez
Department of Mathematics and Institute of Bioinformatics, University of Georgia, Athens, GA, 30602, USA.

Abstract

Muscular dystrophy (MD) 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. These dynamics have begun to be elucidated in the last two decades. This article reviews mathematical models of MD and models that could be used to study molecular and cellular components implicated in MD progression. A biological background for these processes is also presented. Molecular effectors that contribute to MD include mitochondrial bioenergetics and genetic factors; both drive cellular metabolism, communication and signaling. These molecular events leave cells vulnerable to mechanical stress which can activate an immunological cascade that weakens cells and surrounding tissues. This review article lays the foundation for a systems biology approach to study MD progression.

Introduction

Muscular dystrophy (MD) describes generalized progressive muscular weakness due to the wasting of muscle fibers. While this is an umbrella term used to describe a wide range of muscle wasting diseases, the two types that have been most extensively mathematically modeled are Duchenne’s (DMD) and Becker’s (BMD) 1,2, with DMD being the most common childhood form of MD in males. Children affected by this X-linked recessive disorder have an average life expectancy in the mid-twenties, typically becoming fully wheelchair dependent by their teens3. Most forms of MD leave striated muscle cells with reduced contractile abilities, leading to a wide range of phenotypic expression in patients from fatigue to drooping eyelids.

Previous research attributes pathogenesis of DMD and BMD to either absent or partial forms of the dystrophin protein4. From a cellular perspective, the basic contractile unit of the muscle is the sarcomere. The dystrophin protein is located between the sarcolemma, the outer membrane of the sarcomere, and outer layer of myofilaments, providing a scaffold for muscular contraction. Weakness typically begins in extremity muscles, propagating in a proximal–distal direction, until ultimately affecting the diaphragmatic muscles responsible for breathing3,5. Hypertrophy of cardiac muscle cells is an additional complication associated with both DMD and BMD. The loss of function associated with both typically leads to premature death6.

Epidemiological impacts of MD remain difficult to pinpoint definitively. Each subset of MD contains its own range and pattern of pathogenesis and progression. We can further classify subsets of MD in (ordered in decreasing rates of prevalence): dystrophinopathies, laminopathies, dystroglycanopathies, sarcoglycanopathies, and alternative congenital forms7,8. In the United States, dystrophinopathies have a prevalence of about 15 cases per 100,000 people9. The following list describes the rates of prevalence of MD per 100,000 cases in Northern England8:

Myotonic muscular dystrophy is the most common adult onset MD. The prevalence rate is about 10.6.

With respect to dystrophinopathies that manifest as dystrophin absence or deficiency, DMD and BMD have been most extensively studied. Incidence of dystrophinopathies is about 8.5.

Facioscapulohumeral muscular dystrophy (FSHD) has a prevalence rate of about 4.

Laminopathies, dystroglycanopathies,and sarcoglycanopathies are typically classified as Limb-Girdle MD (LGMD). The prevalence rate is about 2.3.

Collagen VI deficiencies that take form in Ullrich Congenital MD (UCMD, also known as Ullrich Scleroatonic MD) have a prevalence rate around 0.13.

Most alternative congenital forms have rates below one.

The goal of this paper is to offer a survey of mathematical models relating to MD, as well as identifying potential opportunities for quantitative research that have yet to be explored within the context of MD. In the context of this paper, mathematical models refer to quantitative approaches that start from mechanistic principles and result in simplified representations of the phenomena that is studied. The modeling formalisms include algebraic expressions, ordinary and partial differential equations, discrete dynamical systems, agent-based models, etc. Mathematical methods differ from statistical methods in that the latter are applied posteriori to broad families of data representations, usually to test hypothesis of significance whereas mathematical models can be used a priori. Clearly, we exclude from this paper statistical analysis of experimental data. However, we include in the scope of the present work strictly computational methods such as network reconstruction.

Each section describes either progression or pathogenesis of several forms of MD. Fig 1 illustrates how the sections relate in some forms of MD. At their roots, most forms of MD have genetic causes; the section Genetic Models presents models of gene-gene interactions intended to discover pathways either specific to a given form of MD or shared between multiple forms of MD. Vulnerabilities caused by genetic mutations result in an increased probability of muscle damage due to normal stress-strain of muscle usage; the section Muscle Models considers simulations of stress-strain during muscle contraction. Muscle contraction requires ATP originated in the mitochondria of which dysfunctions have been implicated in MDs; the section Mitochondrial Models discusses models of cell death initiated by mitochondria as a result of muscle cell damage. Both muscle damage and cell death initiate an acute immune response; this response attempts to regenerate the damaged or dead muscle tissue. Apoptotic signals from immune cells during muscle repair further antagonizes the mitochondrial response. In a patient with some forms of MD, the mitochondrial and immunological responses cause extended periods of inflammation; chronic immune responses have been implicated in the deterioration of muscle fibers in MD patients. The section Immunological Models introduces models of inflammation — both acute and chronic — and tissue repair.

Fig. 1: Conceptual linkage between the sections of this article.

Most types of MD are genetic disorders. Vulnerabilities due to genetic mutation result in increased stress-strain on muscles during contraction. Mitochondrial and immunological responses continue for extended periods; this causes a chronic immunological respone that perpetuates MD progression.

Genetic Models

Most types of MD are believed to have genetic origin4. The genes implicated in MD have functions in numerous biological processes, thus obscuring a univocal characterization of MD pathogenesis. Although each form of MD is caused by different genes, MDs often share similar phenotypes. Finding genetic pathways shared by different MDs could allow for interventions at the downstream node. Furthermore, studying these pathways allow for a greater understanding of how these genes normally interact.

The genetic causes of observed phenomena can be inferred computationally through the reconstruction of gene regulatory networks (GRNs). The characterization of a GRN usually starts with experimental measures of many genes at at least two conditions or points in time. Gene regulatory networks (GRNs) can provide insight regarding subtleties in disease development; graph theoretic models like Boolean Networks10, Bayesian graphs11, and Petri Net12 have emerged as important tools for understanding and reconstructing GRNs. To develop models of GRNs, researchers may employ database like OMIM13 or computer programs like Snoopy14 that use experimental genetic expression data to reconstruct such networks. GRNs can then be studied to find nodes that have more parent genes, are up-regulated or down-regulated by parent genes, affect more downstream genes, or express more often by different but related diseases like MDs. For time dependent data or complex processes, though, directed graphs may be impossible to construct. Instead continuous models like dynamic Bayesian networks15 and temporal Bayesian classifiers16 can be constructed. Although continuous models are unable to answer qualitative questions about the GRN, continuous models can estimate quantitative values. Lack of data hampers continuous models ability to explain large complex GRNs.

Published Genetic Models About MD

Relatively few GRNs have been constructed for MD; among the few successful studies are Grunwald et al. (2008)12 and Tucker et al. (2006)16.

Grunwald et al. (2008)12 created a Petri net of the downstream gene affected by dystrophin (DMD and BMD). The model incorporated sixty-four places and eighty-eight transitions. The authors used Muaritius maps to analysis this model to elucidate the important nodes through knockout experiments. The knockout of the RAP2B-calcineurin pathway caused a 78% disturbance of the whole model. This — the authors concluded — made that pathway important for future study in DMD/BMD therapeutics.

Tucker et al. (2006)16 constructed a continuous model of a GRN involving DMD, LGMD2C, and LGMD2E by using temporal Bayesian classifiers. The model was trained using mouse models (mdx, bsg, gsg) of each form of MD. the model predicts that genes Dlk1, Dusp13, and Casq2 play a part in all three MDs and should be studied further as a possible influence of MD pathogenesis. The role of Dlk1 and Casq2 as shared pathways in MD has been independently verified by statistical studies of gene expressions by Vinciotti et al. (2006)17.

Modeling Other Implicated Genetic Pathways

The task of elucidating gene regulatory networks is a vast area of research (see e.g. DREAM challenge 200818); given the correct data most of these model schemes could be used to study MD. A modern review on this subject can be found at Chai et al. (2014)19. Fig. 2 illustrates the extracellular and intracellular domains of a muscle cell, relating proteins embedded within and associated with the nucleus, the dystrophin-glycoprotein complex (DGC) and the extracellular matrix (ECM). The nuclear outer membrane (OM) is an integral part of the rough ER (Fig 2). Destabilized lamina release proteins that diffuse along the OM to the ER. To aid in future modeling, Fig. 3 includes the known single nucleotide polymorphisms (SNPs) that have been strongly implicated in MD and the genes affected. Publicly available databases can be used to investigate each SNP’s related genes, proteins, and pathways (e.g. SNPedia20, OMIM13, and 1000 Genome Project21).

Fig. 2: Location of Genes Implicated in MD

A diagram of cell structure proteins implicated in MD pathogenesis. The corresponding SNPs for each type of MD shown here are given in Fig. 3. The proteins affected by many types of MD have yet to be elucidated.

Fig. 3: SNPs Implicated in MD

SNPs Implicated in different types of MD. SNP IDs correspond to SNPedia.com.20

Muscular Models

Although vulnerabilities caused by mutated genes, poorly regulated immune responses, and damaged mitochondria perpetuate MD, muscle wasting is the main phenotypical effect of MD5. Damage to muscle fibers caused by contractions exposes muscle tissue to the effects of chronic immune activation and mitochondrial dysfunction. Fibrosis due to chronic immune activation reduces vulnerable muscle cell’s contractual force and flexibility. Understanding these interactions by using mathematical models could offer a macroscopic view of MD progression.

One muscle cell contains thousands of sarcomeres interlaced and primed for contraction. A single sarcomere is made up of parallel thin and thick filaments called actin and myosin; these entwined filaments pull on each other during contraction22,23 . MD leaves muscles vulnerable to unhealthy contraction which exacerbates damage in a cycle that proves catastrophic to muscle cells perpetuating MD progression24,25.

Fig. 4: Sarcomere Regional Classification.

Muscle contraction occurs when calcium ions allow crossbridge sites on myosin to bind with actin and slide actin into the H zone. This action called a powerstoke is followed by ATP binding to the crossbridge allowing the crossbridge to release the actin; Rebinding and powerstroke may than occur again until calcium ions depart the muscle22,23.

Published Muscular Models About MD

Fibrotic tissue accumulation caused by chronic inflammation, contractures, and vulnerabilities in cell structure due to genetic mutation, complicate the process of modeling muscle tissue in MD. Of those attempted, Virgilio et al. (2015)26 and Zöllner (2015)27 created models describing the macro-scale complexities of fibrotic tissue and contractures. Lammerding et al. (2004)28 considers cell deformation in the cytoskeleton due to lamina-deficiency.

Model of fibrosis: Virgilio et al. (2015)26 created a simulation to study the effects of various pathologies including MD. The simulation used an agent based model to create fascicle geometry with the addition of fibrosis and fatty tissue infiltration followed by the use of the micromechanical model described by Sharafi and Blemker (2010)29. This model, using only in silico methods, agreed with the results that fibrosis aggravates the symptoms of DMD.

Model of contractures: MD patients can experience higher-degree of ankle plantarflexion due to contractures; “toe-walking” can make plantar flexor muscles such as the gastrocnemius, soleus and peroneus therapeutic targets. Although the model concerns with long term use of high heels, Zöllner’s model (2015)27 could be used to describe toe-walking in MD. Due to removal of sarcomeres, fascicles shorten with long term toe walking causing less than optimal motion range and increased muscle stress. This increased stress coupled with DCG vulnerabilities accelerates fibrosis and adipose replacement.

Model of cytoskeleton: Nuclear lamina-deficiencies are linked with cytoskeletal defects in muscle cells30,31. Among the mechanisms for these deficiencies are: (i) Destabilized DGC due to cytoskeletal compromise within the plasma membrane leads to the aggregation of mechanical stress in contracting muscle cells28. (ii) Histology in lamina-deficient diseases is commonly concomitant with maligned and internalized nuclei in muscle fibers30,31. (iii) Chronically strained A/C deficient lamins result in cells experiencing higher rates of apoptosis28,32. Intriguingly, regenerating muscle fibers also share high levels of internalized nuclei- the distinction between pathology and regeneration is unclear presently3132.

Lammerding et al. (2004)28 modeled nuclear misalignment and cytoskeletal stress in lamin A/C-deficient mouse embryo fibroblasts. A/C lamins form compartments for splicing factors as well as RNA polymerase II transcription. Inhibition of A/C laminis suppresses RNA polymerase II-dependent transcription in mammalian cells, while its as a scaffold in nuclear compartments remain unknown. A sinusoidal force with amplitude 0.6 nanonewtons (nN) and with frequency 1 Hz, offset 0.6 nN, was applied through a magnetic trap. Cylindrical coordinates (r,θ) were used to measure bead displacement with the magnetic bead at the origin and θ=0 for the force direction. The equation represents the induced strain field described by the analytic cell mechanics model proposed by Bausch et al. (1998)33 ur is the radial component of the induced bead displacement as a function of the applied force, F, cell stiffness μ* the characteristic cut of radius κ-1, the distance from the magnetic bead center r, and the polar angle θ. Thus,

Fig. 5: Lammerding et al. Model

is the radial component of the induced bead displacement as a function of the applied force, F, cell stiffness the characteristic cut of radius , the distance from the magnetic bead center r, and the polar angle . Assuming the cytoskeleton to be incompressible, we can set .

K0 and K1 are modified Bessel functions of respective order 0 and 1 with:

Fig. 6: Bessel Function of Order 1

Letting for incompressibility and metallic bead contact radius of , parameters and were computed by fitting Fig. 6 to the bead displacement.

In this model, the cytoskeleton was exposed to the same biaxial strain as the membrane; for each cell type, nuclear deformation increased approximately linearly with applied membrane strain. Both lamin A/C deficient cells result in decreased nuclear stiffness and altered nuclear mechanics. Under resting conditions, they found that the integrity of the nuclear envelope was maintained in A/C deficient cells. However, under pressure, the control cells maintained nuclei integrity far more than lamin-deficient cells, though both could be ruptured. Increased vulnerability of A/C deficient cells to mechanical stress comparatively was also concluded, with an increase in both necrotic and apoptic cell fraction28. Necrosis mediated through nuclear rupture is not wholly attributed to vulnerability to mechanical stimulation in this model; only about 3-5% of A/C deficient nuclei ruptured.

Modeling Other Implicated Muscular Pathways

Future muscle models into MD could follow several pathways including crossbridge modeling and whole tissue stress-strain simulation. Models of crossbridges combine the physiology of muscle contraction with ATP production in mitochondria. Stress-strain simulation allows for quantification of damage to tissue that activates the immune response.

Cross-bridge models: The crossbridge cycle allows for the ratcheting of muscles during contraction. To understand crossbridges, we must describe the regions of the sarcomere. The A band is made up of thick myosin filaments; throughout contraction its length remains constant centrally localized to the H zone. Actin filaments are laced with the protein tropomyosin. At rest, tropomyosin covers the hinges that catch on actin. Alternatively, the I band is composed of thin actin filaments that alter its length between myosin filament pairs. Structures called Z discs fetter actin filaments at their opposite ends. In S1 regions, myosin edges are hinged and highly flexible, catching on actin and releasing in an indefinite binding cycle (Fig. 4). Upon their release, myosin filaments perform a power stroke catalyzed by the hydrolysis of ATP. ATP is the means for crossbrige formation; it is the release of phosphate during ATP hydrolysis that contracts the S1 region22,23. Calcium provides the means for binding whereas ATP drives contraction. Upon the release of calcium, the protein troponin pushes the tropomyosin to expose binding sites for actin. Upon exposure and a threshold level of ATP expression, the contraction cycle begins34.

Smith et al. (2013)35 described this procedure using a directed graph. The model uses twenty-five components in four domains to describe (mostly protein) interaction during a muscle contraction. The model is a micro scale representation of muscle contraction since it only incorporates one crossbridge. This model could be used to show how lack of ATP due to mitochondrial defects in MD could effect muscle contraction.

Huxley (1957)36 proposed an early mathematical description of this power stroke process. Distance from binding site to crossbridge is taken as the independent variable while function is the probability of a crossbridge being bound at position at time . This yields a conservation law of:

Fig. 7: Crossbridge Conservation Law

where is the velocity of actin to myosin, is the binding rate of crossbridges at position , and is the unbinding rate of crossbridges at position .

The rate of energy release, , by ATP is:

Fig. 8: Rate of Energy Release by ATP

where is the total number of crossbridges at , and is the energy released by a single crossbridge.

This model captures the key physiological properties of crossbridges but fails to illuminate the biochemical interactions; Pate (1989)37 proposed a model based on the Huxley model rectifying the biochemical interaction issue.

The weakening of cardiac muscles due to fibrotic/adipose tissue replacement and metabolic dysfunction is known as cardiomyopathy. This weakening can cause fatigue, swelling in extremities, arrhythmias, chest pain, and possibly heart failure. Cardiomyopathy is the end result of most moderate to severe forms of MD38,39. Understanding how MD results in cardiomyopathy could save lives.

Tewari et al. (2016)40 expanded on Huxley’s concepts with a model of cardiac muscle that connects metabolic and mechanistic crossbridge actions. Using five state variables, the model allows for experimenting with lack of ATP and Ca2+ in cardiac muscle contraction. Additionally, passive muscle action is modeled using spring dynamics. The model was able to fit several experimental data sets of cardiac contraction under different ATP and inorganic phosphate concentrations. The model does not account for three dimensional dynamics nor fibrotic/adipose tissue influences. An expansion of this model to include calcium has also been created41.

Stress-strain models: Levels of healthy isometric force production can be maintained under the strain of intense activity. Rest commenced with the buffering of calcium levels in the cytosol reduces mechanical stress to a point; inappropriate high or low levels of activity leave muscles vulnerable to damage. Furthermore, both types of activity are associated with increased levels of cytosolic calcium additionally linked with fiber damage and apoptosis42.

A.V. Hill (1938) 43 created a Force-Velocity relationship model to understand isotonic and isometric muscle contractions; his work provides a basis for multibody, dynamic musculoskeletal modeling and simulation.

Fig. 9: Hill Model

relates rate of muscle contraction (shortening length), v, to the load p where and are constants given by experimental data and is the isometric force.

Although Jewell and Wilkie (1958)44 demonstrated that this model lost accuracy when exposed to sudden changes to muscle length, Hill’s model remains pivotal and is integrated in many other models.

Van der Linden et al. (1998)45,46 attempted whole tissue interaction simulations examining aponeurosis/muscle under stress, limited by computational power, they were restricted to 2D and simplified 3D models. Johansson et al. (2000)47 improved upon van Leeuwen-Kier’s (1997)48 model of squid tentacles by separately modeling active muscle attributes with Hill’s force-velocity model and describing passive elements as a hyperelastic material. Unfortunately since they used parameter values based on van Leeuwen-Kier’s research and ANSYS, Johansson et al. failed to significantly improve upon the predictions made by van Leeuwen-Kier. Yucesoy et al. (2002)49 introduced a similar model expanding on Van der Linden’s work. Like Johansson et al., Yucesoy et al. separated muscle into the extracellular matrix (passive) and myofiber (active); from that they created a linked fiber-matrix mesh that fuses a passive element and an active element. This “two domain” approach allowed a glimpse into the interaction of muscle fibers and the extracellular matrix. In an effort to understand the effects of geometries on muscle tissue, Sharafi and Blemker (2010)29 devised a model where actual rabbit muscle biopsies were estimated with linear functions. They captured both the geometries of fibers (microscopic level) as well as the fascicles (macroscopic). Since they were modeling muscle stress, Sharafi-Blemker only modeled the passive elements of muscle which were described as a hyperelastic, nearly incompressible material. In opposition to assumptions made by earlier modelists, Sharafi-Blemker discovered that fascicle display anisotropic characteristics.

Mitochondrial Models

Mitochondria are organelles inside a typical cell. Mitochondria produce ATP, which is the currency for sarcomere contraction. In a single cell, there can be hundreds of mitochondria; these organelles are responsible for far more than ATP synthesis. Mitochondria also regulate cell life-death cycle with implications to the immune system. Mitochondria synthesize protein encoded in mitochondrial DNA (mtDNA) in addition to dividing independently as needed inside the cell. mtDNA is vulnerable to mutations that perpetuate the mitochondrial cascade, with rates up to ten times higher than nuclear DNA mutation50,51.

Poorly regulated mitochondria have been implicated in several forms of MD. The distinct role of dystrophin in DMD and BMD’s pathogenesis and progression remains unknown; despite devastating consequences in its absence, the dystrophin protein comprises a mere 0.002% of healthy skeletal muscle proteins52 . Mitochondrial Myopathies are a group of disorders closely related to — and in some articles/websites included as a form of — MD. An accumulation of defective mtDNA could be the cause of these myopathies53.

Published Mitochondrial Models About MD

Mitochondria’s shape, number, and energy processing power are constantly in flux; mitochondrial dynamics alternate between fusion or fission. This fluctuation allows for the dispersal of mtDNA mutation. Tam et al. (2013)54 and Taylor et a. (2003)55 created stochastic models emulating these dynamics.

Tam et al. (2013)54 produced a stochastic model examining rates of fusion-fission that would optimize mitochondrial function and minimize clonal expansion in neighboring mitochondria. This model classifies fusion-fission events in mitochondria using the following criterion:

Fusion events feature nucleoid exchange between mitochondria; one is emptied out to the other and marked for degradation.

Fission sites appear close to the fusion event, regionally containing original nucleoid distributions from precursor mitochondria.

Healthy fission features low levels of exchange of nucleoids, with only mitochondrial matrix contents being mixed.

Larger and longer mitochondria have higher propensities for fission, and smaller ones are more likely to fuse.

The model predicts that protective nuclear retrograde signaling could rescue the mitochondrial cascade through the promotion of mitochondrial nucleoid replication propensity up to sixteen times the basal rate, increasing stochasticity by neutralizing clonal mutant aggregation. Benefits of nuclear retrograde signaling are limited by rates of fusion-fission. Within this simulation, rate of mitochondrial fusion-fission plays a significant role in clonal expansion. Slow exchanges of mtDNA result in homoplasmy, where interventionist retrograde signaling could compound the issue by increasing rate of nucleoid replication. Higher rates of fusion-fission result in a heteroplasmic steady state; increasing levels of mitochondria in cells mix nucleoids faster. These patterns persist regardless of mitochondrial presence in the cell or their replication parameters. A cytoskeletal, cellular model that considers mitochondrial movement independent of fusion-fission, as well as mitochondrial morphology in a differentiated cell context, is needed to further conclude potential therapeutic benefits of retrograde signaling54.

Table 1: Propensity Equations for Tam et al. (2013)

This table lists the equations used in Tam et al. (2013)54 to update the number of mtDNA in the mitochondria where rmax+1 – is the maximum copy number for amplification of mtDNA, Mmito– is the number of mitochondria in a cell at a given time kD – the rate of autophagy, and VF,max – is the maximum propensity of fission.

Propensity

Propensity equations

Upregulated propensity for nucleoid replication due to higher ratio of mutant mtDNA

Propensity for mitochondrial autophagy

Propensity for mitochondrial fission

Taylor et al. (2003)55 also created a probabilistic model. Using this model, the authors were able to estimate mtDNA mutation rate as approximately 5·10-5 mutations per genome per day. This rate guarenteed that clonal expansion in stem cells after 80 years matched known literature. The model provides a simpler example simulating similar results to Tam et al.

Modeling Other Implicated Mitochondrial Pathways

Future research into the effects of mitochondrial dysfunction could follow several paths including models of apoptosis and ATP productions. As discussed above, a newer model studying the effect of accumulations of deleterious mtDNA could be considered for both Mitochondrial Myopathy and DMD. As chronic inflammation has been implicated in several forms of MD, models that simulate apoptosis and necrosis could help explain how acute inflammation turns chronic. Understanding ATP production would also be useful since all muscle contractions require ATP.

Apoptosis models: Byproducts of mitochondrial metabolism include small amounts of electrons that leak from inner membrane complexes and attach themselves to oxygen, forming free radicals called reactive oxidative species (ROS). In a healthy immune response, free radicals such as superoxide and nitric oxide are produced by macrophages for destruction of foreign species. Small amounts of the free radical superoxide produced by the mitochondrion are neutralized by antioxidant enzymes such as superoxide dismutase. Mutated mtDNA leak more ROS in a degenerative mitochondrial cascade essentially poisoning vulnerable cells through ROS release50,51.

Intrinsic apoptosis results from a stressed cellular response. Severe ROS damage can result in cellular necrosis whereas the release of cytochrome c to the cytosol from the inner membrane of the mitochondria triggers intrinsic apoptosis. Depolarization in MOMP (mitochondrial outer membrane permeabilization) in stressed cells triggers the MOMP cascade. MOMP stress markers in survival-apoptotic dynamics feature members of the Bcl-2 protein family and BH-3 only proteins56,58. When triggered by MOMP, mitochondria release cytochrome c and Smac (second mitochondrial-derived activator of caspases) into the cell’s cytosol. By binding with XIAPs (x-linked inhibitors of apoptosis), Smac allows for cytochrome c, Apaf-1, and ATP to combine into apoptosome and cleave procapase-9 forming the initiator caspase-9; cleavage of procaspase-3 by caspase-9 then forms the executioner caspase-3 which in turn activates executioner caspase-6,7 and creates a positive feedback loop by cleaving more caspase-9. These executioners finalize the death of the cell. Apoptosis will not occur if threshold levels of effector capases are not reached56,57,58.

Extrinsic apoptosis occurs when an extracellular self-destruct order is given. Stressed cells signal macrophages to engage with apoptotic cells through phagocytosis as a protective mechanism. Extracellular death ligands act as messagers of these orders by binding with FAS (CD95) death receptors. Subsequently, FAS receptors cluster allowing the binding with FADD (FAS-associated death domain). Recruited by FADD, multiple procapase-8 compile and mutually cleave forming capase-8. This new protein cleaves pro-apoptotic, BH3-only protein Bid forming tBid (Truncated Bid). These interactions lead up to MOMP activation and subsequent cascade; the MOMP activation pathway bridges extrinsic and intrinsic apoptosis. Multidomain proteins are activated by apoptotic tBid activation, which can be inhibited due to protective Bcl-2 proteins56,58.

Fussenegger et al. (2000)56 proposed a model simulating apoptosis to study caspase activation and inhibition (See Table 2). The model confirms experimental observations that Bcl-2 above a critical level effectively inhibits procaspase-9 activation but fails to adequately inhibit procaspase-8 activation, and suppression of FADD’s binding to FAS/FASL complex blocks caspase-8 activation but has little effect on caspase-9 activation. The model assumes isotropic reactions with a well mixed single domain and omits proteins including Bid/tBid, reactions like caspase-8 cleaving of Bid, and bundles executioner caspases 3,6,7 into a single variable. Furthermore, intrinsic and extrinsic apoptosis were not distinguished. Since Fussenegger et al., several models have been proffered to redress omissions.

Albeck et al. (2008)59,60 introduced a model concentrating on the extrinsic apoptosis death switch as well as MOMP interactions. The model represents both cytosol and mitochondria as two separate domains interacting after MOMP with parameters trained by live-cell imaging of HeLa cells. Similar to Fussenegger et al., Albeck et al. bundles many proteins with similar properties — such as caspase-8 and -10 are represented as a single variable C8 — to simplify the model and all reactions are isotropic. However, unlike the previous model, Albeck et al. incorporates a time delay mechanism to compensate for the delay of death ligand reception to MOMP as oppose to the quick death of the cell post-MOMP. Intriguingly, western blot fails to show enough XIAP pre-MOMP to properly inhibit caspase-3. Albeck et al. concluded that another protein/reaction must exist to account for this discrepancy. The model did confirm that MOMP occurs after proapoptotic Bcl-2 proteins reach a certain level depended on the physiological state of the cell. An alternate stable state — partial cell death — is predicted by the model.

To eliminate the isotropic assumption, Huber et al. (2010)57 combined previous models with one-dimensional diffusion PDEs. The typical mass reactions and kinetics are extended by a PDE (Table 2) where vn(x,t) is the chemical reactions. Although their goal was to investigate anisotropic reactions in MOMP, the reaction-diffusion equations remain applicable to wider studies.

Metabolism model: Smith et al. (2013)35 created a directed graph model for ATP production. The model describes the interactions of fifty-nine genes, signals, and transcription factors in 4 domains. Indication of activation or inhibition is included in each interaction although no time scale is given. This network lays the groundwork for future modeling using ODEs and PDEs.

Immunological Models

Immunology describes the coordination of innate and adaptive immune systems to redress threats of the pathogenic, viral, and parasitic kind in addition to the filtration and disposal of excessively damaged cells. Unfortunately, while lifesaving and essential for protection, immune system responses can include self-induced attacks in areas of the body that become so damaged they leave immune cells unable to recognize their healthy peers1,61.

The innate immune system has been implicated in perpetuating multiple forms of MD62,63. Vulnerabilities from missing proteins that stabilize the DCG promote a switch from acute inflammatory response to chronic response; the results of extended inflammation are adipose tissue replacement and muscle fibrosis1,64. This has been confirmed with microscopic analysis of damaged muscle cells lacking dystrophin. Cellular analysis further indicates abnormal, elevated, levels of macrophages, helper and cytotoxic T cells in damaged tissues61.

Published Immunological Models About MD

Cyclic immunological activation purported by mechanical stress is associated with coexisting restorative and degenerative processes in muscle cells; immunological players act both as starters and finishers of apoptosis. CD8+ (cytotoxic T-cells) initiate apoptosis in compromised cells, and the cellular remains of apoptosis are eaten by macrophages61,65. Two models, Dell’Acqua and Castiglione (2009)1 and Jarrah et al. (2014)2, have been proposed using a predator-prey system to study chronic immune activation by employing a log-normal distribution:

Fig. 10: Lognormal Equation

This equation is used to simulate damage to muscle tissue in the following models where is the time to reach peak damage, is the standard deviation of damage, and the magnitude of peak damage at time .

Dell’Acqua and Castiglione (2009)1 is generated by five ODEs in addition to a conservation law (Table 3) describing the immune response of DMD in the mice model. They used COPASI’s66 optimization methods on experimental mdx mice data from Spencer et al. (2001)61 and Wehling et al. (2001)62 to find the best fit parameters. The set of equations used is

Table 3: Dell’Acqua-Castiglione’s Equations

where , , and are kinetic parameters. is a lognormal distribution describing damage to normal tissue. The final equation provides a conservation of tissue.

Equations

Variables

M = concentration of macrophages

H = concentration of CD4+

C = concentration of CD8+

N = percentage of normal muscle fibers

D = percentage of damaged muscle fibers

R = percentage of regenerating muscle fibers

Jarrah et al. (2014)2 (Table 4) refines Dell’Acqua-Castiglione’s model by adding an additional ODE which allows the conservation law to be implicit. The parameters of this model were derived from recent experimental data of mdx mice. Both models assume that the missing dystrophin in the muscle causes damaged muscle cells to initiate the immune response which contributes to their own damage until eventual apoptosis. The set of equations used is

Table 4: Jarrah et al. Equations

where , , and are kinetic parameters. is a lognormal distribution describing damage to normal tissue. This system is similar to the above except for the lack of an explicit conservation of tissue equation. However, the system still has an implicit conservation of tissue with .

Equations

Variables

M = concentration of macrophages

H = concentration of CD4+

C = concentration of CD8+

N = percentage of normal muscle fibers

D = percentage of damaged muscle fibers

R = percentage of regenerating muscle fibers

and the definition of variables is the same as Table 3. Initial conditions have levels of cytotoxic T cell levels at 0; this changes when the impulse damage represented by equation Fig. 10 sets the system into motion. When h=0, the impulse damage is negated and the system remains in a stable state; allowing h>0, the model ensures that T helper cells draw cytotoxic T cells to the damaged region. Damage caused by the immune system reaches a peak in weeks four through eight until the presence of the players wanes. By week twelve, the decreased presence of macrophages, CD4+ and CD8+ T cells (week fourteen) results in diminished levels of degeneration and restoration.

Both models display regions of bistability. Depending on the initial damaged caused and M0, the system collapses to healthy muscle stability or approaches a stability with heterogeneous mixtures of healthy and damaged muscle. This suggests that immune response to muscle damage could be a major contributor to DMD’s pathophysiology1,2 which has been shown experimentally61,62.

These models indicate that the strength of the immune response and maintenance of the positive feedback system relies upon moving past these threshold points to enter another stability state. Driving the system into these recovery regions could prove to be a potential therapeutic target for redressing the role of the inflammatory response.

Modeling Other Implicated Immunological Pathways

The above models incorporate only cellular interactions in an chronic immune response; any multi-scale model would need to include molecular signals that are used in an acute and chronic immune responses to muscle damage. Veltman et al. (2017)67 and López-Caamal et al. (2012)68 introduced models describing molecular signals.

Inflammasome Model: Inflammasome has been implicated in initiating inflammation after an injury or infection occurs by stimulating caspase-1 and activating IL-1β. Both LGMD2B and DMD have shown upregulation of inflammasome69.

A recently published model, Veltman et al. (2017)67, describes inflammasome interactions by using an ODE model to study molecular signaling input to IL-1β activation output. The model also includes inhibitors of inflammasome — IL-10 (signaled by IFNβ) and IFNγ. The model illustrates that negative feedback by IFNγ has less effect on IL-1β than the inhibition by IL-10/IFNβ. This result matches experimental research. The model does not account for anisotropic dispersion nor cellular interactions.

AKT/mTOR Model: Another important avenue for research would be propagation of immune cell signals. López-Caamal et al. (2012)68 created a model describing AKT/mTOR pathway activation by IGF-1. IGF-1 is an important protein released by M2 macrophages during muscle repair that promotes muscle growth, and AKT/mTOR has been implicated as the main signal pathway to promoting muscle growth. The model simulates a positive feedback and diffusion of AKT/mTOR signals once activated. The model also is used to speculate on consequences of signals from other sources like AMPK or pharmaceuticals.

Areas of Future Quantitative Research

Characterizing pathology and pathogenesis of MD requires further study; future models could accelerate therapeutic discovery by testing potential pathways in silico as well as detecting new therapeutic pathways. There are areas for computational research that have not yet been explored mathematically. These pathways may be worthwhile exploring to better understand MD disease development and providing treatments that may delay onset or progression.

The creation of rAAV/AAV1 (Recombinant Adeno-associated virus, Adeno-associated virus 1) delivery system and CRISPR/Cas system indicates the possibility of a genetic cure for MD70,71. Heller et al. (2015)72 overexpressed the human α7 integrin Gene, ITGA7, using the AAV1 delivery system in mdx mice. ITGA7 is a skeletal muscle laminin receptor (Fig. 1) whose overexpression does not cause an immune response in mdx mice. Protective benefits of the DGC were restored with ITGA7 overexpression; lifespans were also prolonged. Xu et al. (2015)73 used CRIPSR/Cas9 to remove mutated exon 23 with the dystrophin genomic region to restore dystrophin expression in the DGC of mdx mice. Clinical trials are currently taking place for LGMD2D74. Mendell et al. (2012)75 wrote a review outlining future work in gene therapy. Both systems require a relatively small number of injection sites. Population dispersion models and GRNs could help in the development and effective administration of both rAAV/AAV1 and CRISPR/Cas systems. Population dispersion and diffusion models could be used to predict the outcome from a series of injections to the spread of the corrected genes throughout the body. Potentially, these models could indicate the most effective injection sites. Beyond the correction of defective genes, both systems could target genes that promote muscle growth and regeneration76. GRNs may help in developing new therapies by finding pathways that both system could target and thereby accelerate the body’s muscle regeneration. These types of treatments could be applicable to both MD patients and in sports medicine.

Myostatin — a TGF-β protein — has long been recognized as a possible therapeutic option for MD. Early murine testing for several MD phenotypes produced mixed results. For DMD (mdx mice) and LGMD2F (scgd{-/-} mouse), Parsons et al. (2006)77 concluded positive results for mice treated early in development before widespread necrosis occurred. Treating sgcg{-/-} mice (modeling LGMD2C) resulted in positive muscle physiology including increase fiber size, muscle mass, and grip force in addition to reduce frequency of apoptosis; however, muscle histology remained unfazed signifying lack of pathology change78. A possible solution proposed by Rodino-Klapac (2009)76 used AAV1 to genetically edit Follistatin (the major myostatin inhibitor). A few experiments show that the new treatment has few side effects and shows similar improvements as other myostatin inhibitors in LGMD2A (Calpainopathy). Future research is needed to show if myostatin inhibition is a means to maintain pathophysiology in MD patients. GRNs and physiological muscle models could be used to understand effective usage of myostatin. GRNs could discover new inhibitors and enzyme activators of myostatin; this may allow targeted genetic editing to regulate myostatin similar to Rodino-Klapac. GRNs might also explain why these methods will work with some types of MD but fail with other types. Physiological models may demonstrate the increase in fiber size, muscle mass, and grip force due to myostatin therapies.

A plethora of models could be used to describe cellular, molecular, and pharmaceutical interactions of the immune system. Precise molecular, mathematical models bridging arginine metabolism with oscillations in macrophage phenotypic expression could be used to model nitric oxide mediated cytotoxicity as well as fibrosis during satellite cell proliferation. Mechanical models for macrophage infiltration and molecular models for macrophage phenotype oscillation could also be useful to better characterize chronic immune system activation due to structural defects. Integration into musculoskeletal simulation may be useful to model the immunological role in MD pathogenesis and progression. Agent based models could be used to imitate immune cells.

MD disease progression also results in alterations to pathophysiology such as gait and muscle atrophy. Noninvasive studies with patients, especially children, could be critical to create a staged model for gait devolution and morphology. Quantifying degrees of eccentric contraction using musculoskeletal simulation could possibly explain selective degeneration in DMD79.

Discussion

With new developments in computational power and data availability, a growing amount of research is using a systems biology approach to understand pathogenesis and progression of disease. Effective and integrated in vitro and in silico models could inform biological phenomena, even without the need of a living subject. For instance, over the last few decades, collagen hydrogel with muscle derived cells (CHMDCs) have promised to revolutionize in vitro experiments and tissue engineering. For CHMDCs to reach the envisioned use, verification by use of mathematical simulations are needed. Recently while examining shape and design, Hodgson (2015)80 used a combination finite elements and agent based analysis to illustrate the lines of principle strain and cell migration in CHMDCs confirming earlier in vitro work by Smith et al. (2012)81 . As MD is a rare disorder, the use of mathematical models could help elucidate the underlining mechanisms of the disease that might not be easily detectable given the limited subject pool.

Although genetic studies have implicated genes as the cause of many types of MD, relatively little is know about about common pathways between these genes that may affect pathogenesis and create similar phenotypes; this necessitates the use of mathematical models describing GRNs (Gene Section). Common genes like Dik1, Dusp13, and Casq21623and there downstream pathways provide future prospects for therapeutic intervention.

Mathematical models of muscles under contractile stress are essential to understanding the long term development of most types of MD. Partial differential equations and agent based models of anisotropic strain from contraction display where cellular rupture and immune response will likely occur. Physical therapeutic and pharmaceutical interventions can be targeted to such areas of high stress to stymie MD progression. Expanding cellular models of the immune response and combining with molecular signals could create a more comprehensive view of muscle tissue regeneration and damage caused by chronic inflammation. Future research could incorporate multiple levels of models into a unified simulation to give a whole view of the progression of MD.

Future research that sheds light on MD disease dynamics, which is likely to occur through mathematical modeling, will provide the means to engage and perhaps ultimately bypass biological systems coordinating together to exacerbate degeneration.

Funding

The authors received no specific funding for this work.

Competing Interests

None of the authors have a financial conflict of interest. M.H. was differential diagnosed (and genetically confirmed recently) with Limb-Girdle Muscular Dystrophy in 2004. A.C. was diagnosed with Progressive Mitochondrial Myopathy in 2010 through a muscle biopsy and differential diagnosis. She worked for the Foundation for Mitochondrial Medicine since 2013 as an intern and now is a volunteer as needed.

Data Availability

This review article does not include any data. All information underlying this study is within the paper.

Corresponding Author

Juan B. Gutierrez, jgutierr@uga.edu

Author Contribution

J.G. conceived of the paper and advised A.C. and M.T. on the organization of the manuscript. A.C. conceived of the mitochondrial subsection and created the sarcomere regional classification graphics. M.T. conceived of the muscle subsection and heavily revised the manuscript after reviewer suggestions. Both A.C. and M.T. contributed to the genetic, mitochondrial, immune, muscle, areas of future quantitative research and discussion sections. All authors reviewed the manuscript.