A Mechanistic, Stochastic Model Helps Understand Multiple Sclerosis Course and Pathogenesis

Heritable and nonheritable factors play a role in multiple sclerosis, but their effect size appears too small, explaining relatively little about disease etiology. Assuming that the factors that trigger the onset of the disease are, to some extent, also those that generate its remissions and relapses, we attempted to model the erratic behaviour of the disease course as observed on a dataset containing the time series of relapses and remissions of 70 patients free of disease-modifying therapies. We show that relapses and remissions follow exponential decaying distributions, excluding periodic recurrences and confirming that relapses manifest randomly in time. It is found that a mechanistic model with a random forcing describes in a satisfactory manner the occurrence of relapses and remissions, and the differences in the length of time spent in each one of the two states. This model may describe how interactions between “soft” etiologic factors occasionally reach the disease threshold thanks to comparably small external random perturbations. The model offers a new context to rethink key problems such as “missing heritability” and “hidden environmental structure” in the etiology of complex traits.


Introduction
Multiple sclerosis (MS) is an immune-mediated disease of the central nervous system with a relapsing-remitting course in the majority of the early stages of the disease [1]. As for other multifactorial diseases, there is no comprehensive overview of the events that lead to the disease. This limits the opportunities provided by the advancements in genetics, immunology, and neurobiology since it is difficult to contextualize each single discovery. The uncertainties in the interpretation of genome-wide association studies (GWAS) reflect, to some extent, this problem. These studies carried the expectation to define the heritable component in multifactorial diseases and, through this, also sketch the nonheritable (environmental) contribution to the phenotype. As largely witnessed by the debate about "missing heritability" in multifactorial diseases, also this powerful approach appears to be in need of interpretative keys as neither genes nor the environment seem to harbour factors that, alone or jointly, are strong enough to explain the disease etiology [2,3].
Likewise situations are rather common in the physics of nonlinear systems; here the observed large variations are explained through the effects induced by small random perturbations [4][5][6]. An example is the theory of the Earth's climate: the cooperative effect of a small stochastic perturbation and periodic forcing (variation of astronomical parameters) provides an amplification of the climate response known as "stochastic resonance" that leads to the transition from a temperate climate to an ice-covered Earth state and vice versa [5][6][7]. 2 International Journal of Genomics Random molecular events can contribute to phenotypic diversity in isogenic populations. Such variation can arise from random fluctuations in gene expression [8,9], and in fact the term "gene expression noise" is typically used in broad reference to variations among seemingly identical cells experiencing the same environment [10]. Noise in genetic circuits can have profound implications also in multicellular organisms enabling physiological regulation mechanisms, differentiation strategies, and facilitating adaptation and evolution [11], often with a qualitatively different outcome than a deterministic one [12]. As such, these discrete and random fluctuations in the expression of individual genes may also participate in the transition between health and disease by exposing hidden genetic and environmental risk.
We show that, in MS, transitions between remissions and relapses (i.e., states of relative health and disease, resp.) indeed occur randomly suggesting that stochastic events may contribute to disease pathogenesis. This is described in a model that accommodates the erratic course of the disease and concedes that "soft" heritable and nonheritable perturbations can be amplified over time. The model is based on the equation for a double-well forced by a random perturbation (noise) [13]. We assume that the patient has two states, disease (relapse) and health (remission), and that the barrier to overcome for going from one state to the other is determined by his genetic background and environmental exposures. The small random perturbation is thought to be the inherent variability in gene expression that triggers the erratic sequence of relapses.
In the likely hypothesis that the mechanisms leading to the transition between remission and relapse recapitulate those that are responsible for the transition between health and disease at onset, time and noise appear as crucial variables that can "amplify" gene-environment interactions and effects until, at one time or the other, the transition towards the disease state occurs in "any" susceptible patient, provided that enough observation time is given.

Patients Data.
We used a dataset of relapses and remissions of 70 patients (28 males and 42 females) with definite MS [14,15] who were monitored prospectively along the years at the MS Clinic of Sapienza University of Rome (Italy) by four experienced MS neurologists. All patients were free of any disease-modifying therapy (all data are antecedent 1993) and had received short courses of corticosteroids if deemed necessary during relapses. Times to disability endpoints were not different between this cohort and published data on the natural history of the disease. The period of interest starts with the first relapse at onset and ends with the last one before the shift to a secondary progressive form. A relapse was defined [14] as the occurrence of new symptoms, the reappearance of former ones, or the worsening of current symptoms of at least 24 hours duration but less than 6 month. Its duration was calculated as the interval between onset of the first symptom or sign and maximum improvement of the last symptom or sign. As part of the standard procedures of specialized MS Clinics, which require to firmly establish the presence or absence of disease activity and progression in order to ensure the most appropriate therapeutic interventions, symptoms were routinely verified during visits at the outpatient service of the MS Clinic. Being this a analysis of anonymous clinical data (collected until 1993, before the institution of ethics committee in Italy in 1998), stored for both clinical and research purposes in the database of the University hospital, ethics committee approval and written informed consent of patients are not needed (http://www.garanteprivacy.it/garante/doc.jsp?ID=1884019).
Each subject has been represented by a sequence of plus one (+1) and minus one (-1) corresponding to the states of relapses and remissions, respectively. For clarity, hereafter these two states are denoted as no health and health. Events are reported on a weekly scale (shorter exacerbations have been rounded up to one week).

The Model.
First, the dataset is analyzed to study the time evolution of relapses, their time duration, and their statistical distribution. Secondly, a simple model consisting of a mechanistic component and a stochastic perturbation [13] is introduced to represent the main characteristics observed in the clinical histories of patients. The mechanistic component provides the deterministic (predictable) component of the model, while the stochastic perturbation describes an "external stimulus" whose behaviour is intrinsically nondeterministic (random) and that can be analyzed only in probabilistic terms. The latter is usually called "stochastic forcing" since it acts as a force that perturbs the deterministic component, making the model solution unpredictable.
The underlying hypothesis of the model here proposed is that the health state of a patient in time behaves likewise the motion of a particle of unit mass in a double-well energy potential as schematically illustrated in Figure 1. Now, given any initial position of the particle and if no perturbation is applied on it (i.e., under deterministic conditions), the particle will fall in one of the two potential wells and will remain there for infinitely long time, independent of whether the well is deep or flat (Figure 1(a)). If, instead, there is a strong stochastic perturbation acting on the particle, it will randomly move forth and back between the two wells, independently of their depth. The most interesting situation for us occurs when there is a small/moderate stochastic perturbation acting on the particle: the particle will typically stay for some time in the well it occupies, until the random diffusion drives it over the potential barrier into the other well ( Figure 1(b)). It is intuitively clear that the exit from a flat well happens faster than the exit from a deep well.
In the phenomena of diffusion, like the one just described, the random perturbation is represented by the Wiener process [16], which consists of a normally distributed white noise (i.e., the amplitude of the noise is normally distributed with given mean and standard deviation, and all frequencies are involved with no time periodicity).
We note that the noise in physical systems was considered for a long time just as a nuisance that degrades the signal. Later on it was found that, in complex (nonlinear)  Figure 1: Schematic picture of a particle in a double-well energy potential. (a) The unperturbed particle falls into one of the two wells from its initial position (deterministic case). (b) The particle subjected to a random perturbation jumps from a well to the other randomly (stochastic case).
processes, as the two-state process here considered, it may introduce a longitudinal dependence corresponding to the erratic switching between the steady states that, otherwise, would be persistent forever [5]. This property of noise should be distinguished from the mechanism of "stochastic resonance" introduced to explain the cyclic recurrence of ice ages and that is based on the combined effect of noise and a periodic external forcing [6]. In that case, useful applications have been found in physical, technological, and biomedical contexts ( [17][18][19][20] and references therein).
Coming back to the metaphor particle-patient, the noise may represent small biological variations of random character experienced by the patient at each time. The bottom of the two wells represents the health and no health states, while the height of the wells is the barrier to be overcome for going from one state to the other, here supposed to be determined by the genetic background and environmental exposures of the patient. Only the effect of noise over time can let the patient change his state of health into disease and be subjected to random relapses.

Some Evidence from Clinical Data.
By analyzing the dataset of 70 patients, some evidence emerges. Figure 2 shows the time evolution of the health states for three sample patients. The time series suggest the unpredictable nature of relapses (+1 states) followed by periods of relative quiet remission with no new signs of disease activity (-1 states). Also, the time in a health state seems much longer than the one spent in a disease state. It is worth noting that the transition between the two states, by construction, is here represented by a step function (abrupt change of state). Such a choice can be reductive or unrealistic since patients usually report a slow or subacute onset of the relapses. However, this feature does not affect the results presented below.
As shown in Figure 3, the duration of the relapsingremitting phase (i.e., period of analysis) varies from one patient to the other, from a minimum of 40 weeks to a maximum of 1311 weeks (about 27 years). As illustrated by Figure  Similarly, the distributions of the time duration (or exit time) of relapses and remissions follow exponential decaying behaviours (Figures 4(a) and 4(b)). The mean duration of disease attacks is about 4.3 weeks, with a minimum of 1 week to a maximum of about 24 weeks (rare events). The mean time spent by patients in the health state is instead about 100 weeks, with a minimum time interval between two consecutive disease events of a few weeks up to about 1000 weeks in exceptional cases. It is worth noting that the absence of any peak in the distributions at some specific time corroborates what has been noted before for sample patients ( Figure 2): disease events occur randomly in time. Moreover, the mean duration of disease events is much shorter than that of remission.
In the following, we will show how these features can be taken into account in a simple mechanistic model forced by random perturbations. The aim is to contribute to the understanding of the mechanisms that underlie the disease, highlighting the role of noise in a system (patient) that is primarily governed by mechanistic laws.

The
Model: Mechanistic Component. The model to simulate the random transitions between the health and no health states (-1 and +1, resp., as in clinical data) is here introduced. Let be the health state of a patient (the position of the particle in our metaphor). For a two-state process, as the one we are interested in, the time evolution of is given by the following equation: where is the time, is a control parameter, and and are the infinitesimal variation of the health state and of time, respectively. The above equation has three steady (independent on time) states. They are 0 = 0, 1 = −1/√ , and 2 = +1/√ . Note that in case = 1, we have 1 = −1 and 2 = 1. Now, in agreement with the convention used for clinical data, let 1 be the state of health and 2 the state of no health.
The potential , associated with (1) can be easily computed [6] as, It can be noted that ( ) is a symmetric function with respect to a change of sign of the variable (it is the same for greater than zero or less than zero), describing a double-well. For clarity, we show in Figure 5 the potential as a function of for two values of the control parameter: = 1 and = 0.7. It should be noted that the difference of the potential in 1 ( 2 ) and 0 (denoted by Δ in the figure) is the barrier that must be overcome to jump from one state to the other; it increases when is decreased with respect to the reference value = 1 (the opposite occurs when is greater than 1). In the context of MS, such a barrier is though to be set by the combination of heritable and nonheritable risk.
Thus, according to (1) and (2), we have (a) if, by any chance, at initial time the patient is nearby 1 , he will be there forever; (b) if, by any chance, at initial time the patient is nearby 2 , he will be there forever; International Journal of Genomics   Figure 5: Potential (double-well) as a function of the state variable . The potential associated with the mechanistic model described by (1) is plotted for two values of the parameter controlling the height of the wells: (a) = 1 and (b) = 0.7. The variable represents the health state of a patient (i.e., the position of the particle in our metaphor). The health state is denoted by 1 , the no health state by 2 , and 0 is the unstable steady state between them. Δ is the barrier whose height is controlled by the parameter . Dashed lines denote the upper and lower bounds of the potential barrier.
(c) if, by any chance, at initial time the patient is nearby 0 , he will have 50% chance to be forever nearby 1 or 2 . Now, as shown by clinical data (Figure 2), the patient ( ) stays a longer time in the health state compared to the no health state. This feature suggests that the depth of the two wells differs and, in particular, the one associated to the health state is deeper than the other, so that it is more difficult to exit from the health condition than from the relapse. In our mathematical model, this implies the introduction of an asymmetry in ( ), as where the parameter determines the magnitude and shape of the asymmetry.

6
International Journal of Genomics Figure 6 illustrates an example of asymmetric doublewell ( = 0.08) for = 1 and = 0.7. In Figure 6, the barriers are denoted as Δ 1 and Δ 2 and are the difference in the potential between the states 1 and 0 , and 2 and 0 , respectively. In all the three circumstances above, such barriers separate the states 1 and 2 , preventing the switching between the wells. Now the question is: can some external random stimuli provide enough energy for diffusing across the barriers, regardless the initial conditions?

The External Random Stimuli.
Unless new risk factors are discovered, those identified so far in MS appear unsuitable as deterministic triggers of continuous and random transitions between health and disease, inducing very different disease courses not only in different patients but also in the same patient at different time points. As noise is a source of variability in biological systems, we tried to model it as a driving force. This kind of forcing must be specified through its statistical properties. Let us suppose that on average the noise has zero effect (i.e., its time mean is zero), has a given variance, and it is time decorrelated (i.e., at a given time the value of the stimulus depends only on the previous time): this is the zero order approximation of the statistics of noise (i.e., the Wiener process previously discussed). By introducing such external random perturbation into (1) where 1/2 is the noise of variance . Because of the external random stimuli, this equation has important deviations from its deterministic version described by (1): now the patient can change randomly his state from health to no health and vice versa.
Averaged over time it is impossible that crosses the barrier since the noise, by construction, has zero mean. However, given the statistical property of the stimulus, in due time, we expect that a sequence of small driving occurrences may be equally signed so that may overcome the barrier and switch from one state to the other. In other words, is harvesting energy from the stimuli for crossing the barrier. For the sake of clarity, the relationship among disease characteristics and model parameters are summarized in Table 1.
As we shall see, this model fits nicely the statistics of the exit times of disease events and explains a few paradoxes encountered in understanding the disease course.

Model Solutions.
To better understand the model, let us look at the solutions of (4) for different values of the parameters. Setting = 1, the solution of (4) as a function of time is shown in Figure 7 for the symmetric and asymmetric double-well cases ( = 0 and = 0.08, resp.). According to our metaphor, the patient alternates periods of wellness and disease randomly and, as expected, when a small asymmetry is taken into account, the patient spends more time in the health state than in the no health state. This means that there is a finite probability that for any initial condition nearby 1 (or 2 ) the patient will jump across 0 at a finite time (i.e., time duration in a given state), switching from one state to the other.
As seen from Figure 4, the duration time of each state, , is a random variable. The mean value of in the health state ( 1 ) can be estimated [6] and is given by where is the Nepero number. It is interesting to note that, knowing the mean values of duration times in the health and no health states, it is possible to estimate the ratio between Δ 1 and Δ 2 , say the asymmetry of the double-well. For the dataset discussed in the previous section, we have 1 ≈ 100 weeks and 2 ≈ 4.3 weeks. On the other hand, taking the logarithm of (5) and (6), and performing their ratio, we get Thus, the ratio of the logarithms of duration times in the states 1 and 2 provides an estimate of the ratio between the barriers that must be overcome to move from one state to the other. By considering the duration times 1 and 2 averaged over all the cases of 70 patients it is found Δ 1 /Δ 2 ≈ 3.1. This implies that the well in 1 is about three times deeper than that in 2 . We point out that the ratio between the depths of the wells can be estimated for each patient, provided that the mean duration times of the two states are known from the clinical history. As an example, we have estimated such a ratio for the three sample patients considered in Figure 2. We found At this stage, the potential associated with each patient can be estimated by fixing = 1 and changing so that the two wells are asymmetric in the predetermined ratio. Plots of the three potentials are shown in Figure 8. By assuming that random stimuli of the same variance are applied to the three patients, the shapes of the potentials reflect the differences in the duration times of health and disease states observed in Figure 1. The first potential has a very deep well associated with the health state and a very shallow well associated with the no health state: this configuration lets the patient be affected by a few disease attacks of short duration. Given the same observation period, the patient no. 32 is expected to  Figure 6: Asymmetric potential as a function of the state variable . The asymmetry in the depth of the wells is obtained setting ̸ = 0 in (3). Here the effect of ( = 0.08) on the shape of the potential is illustrated for the same values of the control parameter used in Figure 5: (a) = 1 and (b) = 0.7. Δ 1 and Δ 2 are the barriers between 1 ( 2 ) and 0 , respectively. Dashed lines denote the upper and lower bounds of the potential barriers. Table 1: Disease characteristics and parameters of the mechanistic stochastic model. The table summarizes analogies and links between MS characteristics (left column) and model parameters in (4), in the right column.
Disease characteristics during the relapsing-remitting phase Mechanistic, stochastic model Two states are observed that we referred to as health and no health.
( ): double-well potential (two minima associated with the two states).
The disease develops in patients with a certain degree of heritable and nonheritable risk.
: control parameter that sets the barrier to overcome for going from one state to the other.
It is found that the time spent in the state of health is much longer than in the no health state.
: asymmetry parameter of the potential.
Random variability of biological processes (i.e. randomness in transcription and translation leading to cell-to-cell variations) [8].
Genetic background and environmental factors alone are not enough to trigger the transitions between the two states.
The noise provides the required energy for the random transitions.
have more attacks compared to the previous patient due to the difference in the depth of the wells. Lastly, the patient no. 53 when compared to the patient no. 32 should experience disease events lasting more time because of the deeper well associated with the no health state.

Discussion
The contribution of the heritable or non-heritable factors that drive MS onset and course appear uniformly too small to explain its etiology, variable severity, and erratic manifestations. We present a mechanistic, stochastic model that describes how the effects of the above forces may be randomly amplified over time. Given a time of observation T (e.g., the average lifespan of a person), patients may experience a transition between the two states with a mean time . Thus, the probability to observe a transition requires that T is much longer than . However, since is randomly distributed, by mere chance, the transition may or may not occur within the observation time T. Within this approach, time is conceived not only as a random variable characterizing disease events but also as the "observation time" needed to detect the onset of the disease in "any susceptible" patient. Stochastic noise in gene expression, through its pervasive effects on virtually all biological processes, may be the factor that amplifies and reshapes the deterministic effects of genetic and environmental risk factors. In both cases, the control parameter is set to = 1, and a stochastic perturbation of the same variance is applied ( = 0.13). Under the effect of the random perturbation, the state variable (i.e., the health state of a patient in our metaphor) jumps erratically from −1 to +1 and vice versa (a). When a small asymmetry is introduced in the potential , the variable still changes erratically its state but spends more time around −1 than around +1 (b). The latter case appears suitable for describing the random occurrences of relapses in MS as observed in clinical data. The model is in accord with a variety of data in MS and may be general enough to explain features of other complex traits. First, time is a key variable in this model but also in experimental autoimmunity where there is evidence of a multistep process made of subtle alterationsresulting from quantitative trait loci variations-that may accumulate with time and end up in susceptibility when a threshold for the occurrence of the disease is passed [21,22]. Furthermore, defects of the immune response that originate from predisposing genetic variants are present before any onset of autoimmunity [23]. These relatively recent observations are in accord with classical epidemiological data suggesting that fairly long time intervals are a necessary factor for the development of MS [24,25]. Lastly, time itself remains an essential component of the diagnostic algorithm of MS [26]. It is therefore not surprising that, if time is a crucial factor to understand the nature of the disease in diagnostic terms, it may also be instrumental for the pathogenetic understanding of the process. Second, our model links the degree of genetic variation with the probability that the transition between health and disease states occurs over time. In fact, MS risk alleles influence age of onset [27] and correlation of age of onset in MS relative pairs is proportional to genetic sharing [28]. Interestingly, the model applies also when individuals share genetic protection as witnessed by our observation of an underrepresentation of MS among Italian twins, possibly due to protective factors shared by twins in a Mediterranean area [29]. Furthermore, as expected, a link between genetic variation and transition between the two states over time seems to be maintained during the disease course: with one exception [30], the vast majority of the studies on the natural history of MS has shown that, while the occurrence of relapses is unpredictable, their frequency and severity tend to be conserved as a longer first interattack interval is associated with a better prognosis [31][32][33][34].
While this is the first time that a model combining deterministic factors and stochastic forcing is applied to the interpretation of a multifactorial disease, determinism and stochasticity have been already shown to fit basic biological mechanisms ( [8][9][10][11][12] and references therein). Here, stochasticity is beneficial to "the individual, the colony, or the species" [35] in that it may exert a physiological regulatory role. It may therefore be feasible to seek strategies that try to revert a "pathogenic stochasticity" to its original physiologic role. In the stochastic resonance, introduced to explain Earth's glaciations [5] and sometimes erroneously associated with the random transitions illustrated above, the model solution is driven by the cooperative effect of the stochastic perturbation and the small periodic forcing, leading to a flip from one state to the other with the frequency proper of the external periodic forcing; the small stochastic perturbation provides the energy necessary for the transition, while the periodic forcing sets the periodicity of the transitions. The stochastic resonance occurs for low frequency of the periodic forcing, say for long time periods. For climate changes that induce glaciation cycles [6], such a period is set to 100000 years that describes the long-term variation of the energy input prescribed by astronomical theory. In the case of relapsing-remitting MS, as clinical data suggest, there is not a typical periodicity of relapses induced by some external factors. However, we can imagine to exploit a small external cyclic component in order to "stabilize" the solution around a given state (the health state) "deamplifying" it with some short periodicity, contrary to what happens in the stochastic resonance, in order to prevent relapses. A nonmutually exclusive approach may be to decrease the variance of the noise (possibly acting on fluctuations of gene expression).
This study offers a new interpretative key of the events that characterize the clinical course of MS. By inference, it also suggests that time and noise amplify and adjust the deterministic influence of genes and environment on the disease. It will not be easy to obtain direct experimental evidence for this model at the disease level. However, basic science data-increasingly supporting the notion that noise provides critical functions at the cell level-suggest that this effort will probably be a rewarding one. Thus far the model may contribute to the understanding of phenotypic variation in health and disease states.

Conclusion
Despite the advent of genomewide association studies (GWAS; an approach to look for associations between hundreds of thousands genetic variations and diseases), many questions remain about the causative mechanisms of multifactorial diseases. In particular, the effect size of the genetic variants identified through these studies explains relatively little of the heritability (the proportion of total variance in a population for a particular measurement, taken at a particular time or age, that is attributable to variation in genetic values) of most complex traits. So far, the same is true for environmental risk factors. Multiple sclerosis is one of the diseases that best exemplify this key problem. By observing the course of this disease in its relapsing-remitting phase we found that a mechanistic model with a random forcing describes the main features of the disease course. If the model is applied to the events that precede the disease (i.e., those that derive from the effects of the genetic and environmental risk factors), it may explain how relatively small effects may be amplified by comparably small random perturbations. As such, the model may inform the design of future studies on gene-environment interactions in multifactorial diseases.