1 Introduction
MS is a chronic autoimmune disorder characterized by immune-mediated destruction of the myelin sheath in the central nervous system (CNS). This demyelination process, together with associated axonal damage, disrupts neural communication and leads to a wide range of neurological impairments. Despite extensive clinical and experimental research, the underlying mechanisms driving disease progression, relapse patterns, and long-term neurological decline remain incompletely understood, particularly in relation to the complex interplay between immune activation and regulatory processes [1].
Mathematical modeling has become an essential tool in biomedical research, providing a systematic framework for analyzing complex biological systems, testing hypotheses, and informing therapeutic strategies. Such approaches have been successfully applied to diseases including cancer, HIV, and infectious diseases, where they have contributed to identifying key regulatory mechanisms and optimizing treatment protocols [2,3]. In contrast, multiple sclerosis has received comparatively limited attention from a mathematical modeling perspective, despite its inherently immune-driven nature and clinical complexity.
This underrepresentation highlights a significant opportunity for the development of quantitative frameworks capable of advancing both theoretical understanding and treatment design.
Existing mathematical models of MS dynamics have primarily focused on specific aspects of the disease, such as T-cell activation, lesion formation, or simplified immune–myelin interactions [4,5]. While these studies have provided valuable insights into individual mechanisms, they generally do not simultaneously incorporate damage-dependent immune activation, regulatory immune responses, treatment-induced modulation, and optimal control within a single analytically tractable framework. The present model addresses these aspects simultaneously, allowing the investigation of treatment-dependent threshold dynamics and long-term immune behavior.
Recent advances in immunology emphasize that MS cannot be fully understood through a purely T-cell–centric perspective, as regulatory immune mechanisms and broader immune interactions play a crucial role in disease progression [6,7]. At the same time, clinical evidence demonstrates that early and sustained treatment can significantly alter disease trajectories. However, there remains a lack of mathematical models capable of quantitatively linking immune dynamics, therapeutic intervention, and long-term neurological outcomes.
Addressing this gap is essential for improving both mechanistic understanding and treatment optimization.
In this work, we develop and analyze a novel mathematical model that integrates damage-driven immune activation, regulatory immune response, and treatment-induced modulation within a unified framework. The model explicitly captures the coupling between myelin degradation and immune activation, leading to a nonlinear feedback structure that governs disease progression. A quasi–steady-state reduction yields a reduced system that remains analytically tractable while preserving the essential biological mechanisms.
The key contribution of this work is the characterization of a treatment-dependent threshold structure that links treatment intensity, immune activation, and long-term disease dynamics. Specifically, we show that the model exhibits two distinct dynamical regimes separated by a critical treatment level: one characterized by persistent oscillatory autoimmune activity and the other by immune suppression and convergence to the immune-free state. Thus, within the proposed modeling framework, treatment governs the qualitative behavior of the system.
More broadly, the purpose of this work is to provide a mechanistic framework for understanding how myelin damage, immune activation, immune regulation, and treatment interact to shape disease progression in multiple sclerosis. Rather than evaluating a specific clinical therapy, the treatment component is incorporated as a generic intervention mechanism that allows us to investigate treatment thresholds, long-term disease dynamics, and the potential impact of treatment on autoimmune activity and myelin preservation.
We further establish fundamental mathematical properties of the model, including existence, uniqueness, nonnegativity, and boundedness of solutions. The reduced system is analyzed to determine the existence and stability of the disease-free, boundary, and interior equilibria. In addition, an optimal control problem is formulated to investigate how time-dependent treatment strategies influence disease burden while accounting for treatment costs.
The resulting optimality system is solved numerically to illustrate the effectiveness of treatment in suppressing inflammatory dynamics and preserving neuronal function.
The remainder of the paper is organized as follows. We first formulate the mathematical model and establish its well-posedness. A reduced system is then derived and analyzed to characterize equilibrium behavior and threshold dynamics.
Next, an optimal control framework is developed and analyzed. Numerical simulations are presented to illustrate the theoretical results and explore the impact of key parameters. Finally, we conclude with a summary of the main findings and directions for future research.
2 Model formulation and assumptions
We begin by making the following biological and modeling assumptions:
The system is divided into three interacting compartments: healthy myelin in the nervous system M(t), autoreactive immune cells (e.g., pathogenic T-cells) I(t), and regulatory immune cells (e.g. Treg cells) R(t).
Autoreactive immune activation is assumed to increase with the level of myelin damage. We define the damage variable by
where Mmax denotes the maximum healthy myelin level and M denotes the current healthy myelin level. Thus, D = 0 corresponds to the absence of myelin damage, whereas larger values of D represent increasing degrees of demyelination. This choice is biologically motivated by the observation that demyelination increases the exposure of myelin antigens to the immune system, thereby enhancing autoreactive immune responses.
Therefore, it is natural to assume that the immune activation rate is an increasing function of the damage level D.
The immune activation rate is modeled by the saturating function
where ηmax is the maximum activation rate and ω is the half-saturation constant. This formulation reflects the biological observation that greater myelin damage exposes more self-antigens, thereby promoting autoreactive immune activation, while the activation rate approaches a finite maximum for large damage levels.
The saturating form of η(D) reflects the biological observation that increasing myelin damage promotes autoreactive immune activation through the release of damage-associated signals, antigen presentation, and inflammatory mediators. However, immune activation cannot increase indefinitely because of finite immune-cell availability, receptor saturation, signaling constraints, and regulatory feedback. Although no direct experimental evidence establishes this exact functional relationship in multiple sclerosis, a Michaelis–Menten/Holling type II response provides a biologically plausible phenomenological approximation while ensuring a bounded activation rate.
Healthy myelin M(t) is naturally regenerated at rate α, and is degraded (i) intrinsically at rate δM and (ii) through immune-mediated attack at rate βIM. Autoreactive immune cells I(t) proliferate at rate η(D)I, are suppressed by regulatory cells at rate μRI, and are removed by natural death and treatment at rate (θ+τ)I, where θ denotes natural cell death and τ reflects therapeutic intervention. Regulatory immune cells R(t) expand proportionally to autoreactive cell presence at rate ρIR and decay naturally at rate λR.
The functional capacity of neurons N(t) deteriorates when myelin falls below its maximal healthy level Mmax. The degradation of neurons is proportional to the damage level D and the current functional capacity N(t).
In the present formulation, N(t) is treated as a downstream measure of cumulative neurological damage. Neuronal repair or recovery, as well as feedback from neuronal damage to the immune dynamics, are not included and are left for future extensions of the model.
Consistent with the biological assumptions stated above, the model is intended to provide a mechanistic description of the interactions among myelin damage, autoreactive immune activation, regulatory immune responses, and treatment intervention in multiple sclerosis. The model structure is motivated by the current biological understanding of disease progression, whereby demyelination promotes immune activation, regulatory immune cells suppress excessive inflammation, and treatment-induced suppression reduces autoreactive immune activity. Although the model does not include all immunological pathways involved in multiple sclerosis, it incorporates the principal mechanisms believed to drive disease progression while remaining analytically tractable.
A complete list of variables and parameters is summarized in Table 1, and the corresponding flowchart of the model is presented in Fig 1.
The model resulting from these assumptions is given by the following system of ordinary differential equations:
where 0<ε≪1 represents the slower adaptation of myelin and neurons compared to the fast immune (I) and regulatory (R) cell populations.
The small parameter ε reflects the separation of biological time scales between the relatively fast immune dynamics and the slower processes of myelin turnover and neuronal decline. This time-scale separation motivates the fast-time reduction used later in the analysis and allows the reduced system to retain the principal biological mechanisms while remaining analytically tractable.
2.1 Ethics Statement
This study is based exclusively on mathematical modeling and numerical simulations. No experiments involving human participants or animals were conducted by the author.
2.2 Use of generative artificial intelligence
ChatGPT (OpenAI) was used to assist with language editing and improving the clarity of the manuscript. The author reviewed and edited the AI-assisted text and takes full responsibility for the accuracy and integrity of the published content. Fig 2 was generated with the assistance of ChatGPT.
2.3 Existence and uniqueness of solution
For the system of ordinary differential equations (1)–(4), we denote the right-hand sides by the vector field
To ensure well-posedness, we assume that F is continuous in t and locally Lipschitz with respect to X on ℝ+4. The presence of the small parameter 0<ε≪1 does not affect these regularity properties. Indeed, all polynomial and bilinear terms (such as IM, RI, and DN) are smooth and therefore locally Lipschitz on ℝ+4.
Moreover, since D=Mmax−M is a linear function of M, it is smooth. The function η(D) is continuously differentiable and bounded for D≥0 (in particular, 0≤η(D)≤ηmax), and hence is locally Lipschitz. Consequently, F(t,X) is locally Lipschitz with respect to X on any biologically relevant domain where 0≤M≤Mmax.
Therefore, by the Picard–Lindelöf (Cauchy–Lipschitz) theorem, for any initial condition
there exists a unique local solution
defined on some interval [0,T1), where T1 > 0 is the maximal time of existence.
2.4 Nonnegativity of Solutions
Theorem 1.Assume the initial data are nonnegative,
Then the corresponding solution satisfies
Proof. We show that the nonnegative orthant is forward invariant.
On the boundary M = 0,
Thus, M(t) cannot cross into negative values.
The I-equation has the form
so
Similarly,
Finally,
so
Therefore, all components remain nonnegative for all t∈[0,T).
2.5 Uniform boundedness of solutions
Theorem 2.All solutions of system(1)–(4)with nonnegative initial data remain uniformly bounded.
Proof. We establish bounds for each state variable.
First, from the myelin equation,
since I≥0. By the comparison principle,
Hence
Next, observe that
satisfies
Therefore,
while
We distinguish two cases.
Case 1:ηmax≤θ+τ.
In this case,
and therefore Let
Then
Hence I(t) is bounded and converges exponentially to zero.
For R(t), we write
Solving this scalar equation gives
Using the bound on I(t), we obtain
Therefore,
Thus R(t) remains uniformly bounded on [0,∞). Consequently, both immune populations are bounded.
Case 2:ηmax>θ+τ.
Consider the Lotka–Volterra comparison system
This system admits the positive equilibrium
Moreover,
is a first integral satisfying
Furthermore,
as
Hence every level set of V is compact in the positive quadrant. Let
and define
Since ΓV0 is compact and the function
is continuous, there exists
Therefore,
where
Since the original immune subsystem is dominated by the above comparison system, its trajectories remain bounded as well.
Finally, for the neuronal capacity variable,
with D≥0. Therefore
which implies
Combining the bounds for M, I, R, and N, we conclude that all solutions remain uniformly bounded for all t≥0.□
Corollary 1.Under the assumptions of Theorems 1 and 2, every solution of system(1)–(4)with nonnegative initial data remains in the region
where
and
Consequently,Ωis a positively invariant and bounded region for system(1)–(4).
3 Fast immune dynamics
For the analysis of the fast immune dynamics, the neuronal variable N is omitted because it is decoupled from the (M,I,R) subsystem and therefore does not affect its equilibrium or stability properties. The variable N(t) is retained in the full model and in the numerical simulations.
The parameter 0<ε≪1 separates the slow myelin dynamics from the faster immune dynamics. We justify freezing M on every fixed fast-time interval 0≤t≤T.
By Corollary 1,
Hence, from
we obtain
Therefore,
uniformly for 0≤t≤T.
Since
it follows that
Moreover,
Thus, using (5),
To compare the full and frozen-myelin immune dynamics, define
and
Then
with
For the same immune state X=(I,R)⊤,
Using I≤Φ* and (6),
where
Since F is locally Lipschitz in X on the bounded invariant region Ω, there exists LT>0 such that
Grönwall’s inequality therefore gives
and consequently
Thus, on every fixed O(1) time interval, the full immune dynamics converge, as ε→0, to the system obtained by freezing the myelin level. Let
The corresponding fast layer problem is
For notational simplicity, the prescribed damage level D0 will henceforth be denoted by D. Therefore, (8) becomes
with D treated as a parameter on the fast time scale.
System (8) is the fast layer problem associated with the full slow-fast system, rather than a conventional quasi-steady-state approximation in which the fast variables are assumed to have reached equilibrium.
In particular, as shown below, when
the positive equilibrium of the fast subsystem is a nonlinear center and is not normally hyperbolic. Thus, the analysis does not require persistence of a normally hyperbolic slow manifold in this regime.
The equilibria and stability properties of (8) are investigated below.
3.1 Equilibria of the fast immune subsystem
The equilibria of system (8) satisfy
and
The immune-free equilibrium is
For a positive interior equilibrium Ei,
and
Hence, the immune subsystem admits the interior equilibrium
provided that
This condition ensures that immune activation exceeds the total removal rate of autoreactive immune cells, allowing the coexistence of autoreactive and regulatory immune cell populations.
At the critical value
the interior equilibrium reduces to the immune-only boundary equilibrium
3.2 Local stability analysis
Let
and
The Jacobian matrix of system (8) is
3.2.1 Immune-free equilibrium.
The immune-free equilibrium of the immune subsystem is
The Jacobian matrix evaluated at E0 is
Hence, the corresponding eigenvalues are
Therefore, the immune-free equilibrium is locally asymptotically stable whenever
and unstable whenever
Motivated by this threshold condition, we define the damage-dependent autoimmune activation number
Using D=Mmax−M, this may also be written as
Consequently,
whereas
The quantity R0 compares the activation of autoreactive immune cells induced by tissue damage with their overall removal rate. It therefore measures the ability of a given level of myelin damage to promote autoimmune activity. When R0<1, immune removal dominates damage-induced activation and autoreactive immune activity declines.
In contrast, when R0>1, damage-induced activation exceeds immune removal, causing autoreactive immune activity to increase. Consequently, R0=1 represents the threshold separating immune control from immune activation.
The threshold condition may also be expressed in terms of the treatment parameter. Defining
we obtain
and
Hence, τc represents the critical treatment level required to suppress autoreactive immune activation for a given level of myelin damage. Since η(D) is an increasing function of D, larger levels of tissue damage require stronger treatment to maintain immune control. This highlights the central role of myelin damage in determining both disease progression and treatment effectiveness.
3.2.2 Interior equilibrium.
Suppose that
so that the interior equilibrium
exists.
Evaluating the Jacobian at Ei gives
The characteristic equation is
Thus,
Therefore, whenever the interior equilibrium exists, its eigenvalues are purely imaginary. Consequently, the linearization predicts a center and linear stability analysis alone is inconclusive.
Nonlinear behavior of the interior equilibrium
Assume that
Define
Then the immune subsystem becomes
The interior equilibrium is
We now show that Ei is a nonlinear center.
Consider the function
Differentiating H along solutions gives
Substituting the system equations yields
Since
and
we obtain
Using
it follows that
Hence,
along every solution trajectory.
Next, compute the first derivatives:
Therefore,
and
Thus, the only critical point of H in the positive quadrant is Ei.
The second derivatives are
Hence, the Hessian matrix is
which is positive definite throughout the positive quadrant.
Therefore, H has a strict local minimum at Ei.
Furthermore,
as
Consequently, every level set
sufficiently close to Ei is a simple closed curve surrounding Ei.
Since H is conserved along trajectories, solutions starting on one of these level curves remain on that curve for all time. Therefore, trajectories sufficiently close to Ei are closed orbits surrounding Ei.
We conclude that the interior equilibrium
is a nonlinear center.
3.2.3 Immune-only boundary equilibrium.
At the critical value
the interior equilibrium reduces to the immune-only boundary equilibrium
Evaluating the Jacobian at Eb yields
whose characteristic equation is
Hence, the point
possesses a double zero eigenvalue, and linearization is therefore inconclusive. This is consistent with the degenerate critical structure described above at
3.3 Global stability of the immune-free equilibrium
Theorem 3.If
then the immune-free equilibrium
of system(8)is globally asymptotically stable inℝ+2.
Proof. The condition R0<1 is equivalent to
Let
Consider the Lyapunov function
Clearly,
for all (I,R)∈ℝ+2, and
Differentiating V along solutions gives
Using system (8), we obtain
Expanding,
Thus,
Since
we get
Moreover.
Therefore, by Lyapunov’s direct method, the immune-free equilibrium E0 = (0,0) is globally asymptotically stable in ℝ+2.
3.4 Treatment-induced stability and long-term dynamics
The threshold quantity
depends on the current myelin level M. Since
and η(D) is an increasing function of the damage level D, it follows that R0(M) is a continuous decreasing function of M. Consequently, there exists a unique critical myelin level Mc satisfying
The slow myelin dynamics provide a biological interpretation of the treatment-induced regimes. When the treatment level exceeds the critical threshold τc, the immune-free equilibrium E0 is globally asymptotically stable and the immune populations converge to zero. In this case, the slow myelin equation reduces to
yielding the disease-free myelin level
When the treatment level is below τc, the interior equilibrium Ei exists and the immune subsystem exhibits sustained oscillatory dynamics around the representative autoreactive burden
Substituting this value into the slow myelin equation gives
which yields the inflammatory myelin level
Since
we obtain
For the parameter regime considered in this work, the disease-free state satisfies
whereas the inflammatory state satisfies
Since R0(M) is continuous and strictly decreasing in M, the unique threshold level Mc must lie between Mi and M0. Therefore,
Defining the corresponding damage levels by
gives
Associated with these damage levels are the treatment values
Since η(D) is increasing in D, it follows that
The threshold τc is the only critical treatment value governing the transition between immune persistence and immune suppression. The quantities τi and τ0 are included only as reference values associated with the representative inflammatory and disease-free myelin states.
Fig 2 summarizes the treatment-dependent stability and threshold structure of the reduced immune subsystem. For treatment levels below the critical threshold τc, the immune-free equilibrium E0 is unstable, while the interior equilibrium Ei exists as a nonlinear center surrounded by closed periodic orbits. At the critical threshold τ=τc (equivalently, ℛ0=1), the system possesses the continuum of boundary equilibria
with
being the limiting position of Ei as τ→τc−. For treatment levels above τc, the interior equilibrium no longer exists and the immune-free equilibrium E0 is globally asymptotically stable. Thus, τc represents the critical treatment threshold separating persistent autoimmune oscillations from immune suppression.
4 Optimal control with smooth treatment profiles
To investigate treatment strategies that reduce autoimmune activity while accounting for treatment intensity and its rate of adjustment, we formulate an optimal control problem over a fixed treatment horizon [0,T].
The constant treatment parameter τ considered in the preceding sections is replaced by a time-dependent treatment intensity u(t). We introduce the control variable a(t)=u˙(t), representing the rate of treatment adjustment. Accordingly, u(t) is treated as an additional state variable, while a(t) is the control to be optimized.
Positive values of a(t) correspond to treatment escalation, whereas negative values represent treatment reduction. Since a∈L∞(0,T) and
the treatment intensity satisfies u∈W1,∞(0,T) and is therefore Lipschitz continuous.
Controlled system and objective functional
The controlled system is given by
subject to
The objective functional is defined by
where
The term AI(t) penalizes autoreactive immune-cell burden, while Ba(t)2 penalizes rapid changes in treatment intensity. The terminal term qII(T) penalizes residual autoimmune activity at the end of the treatment horizon.
The additional linear term Cu(t) penalizes cumulative treatment exposure. Thus, the objective functional balances autoimmune suppression, treatment intensity, and the rate of treatment adjustment.
The admissible control set is defined by
where
The optimal control problem consists of finding a*∈U satisfying
Existence of an optimal control
Theorem 4.Assume that the biological parameters are positive, the initial conditions are nonnegative, and
Then the optimal control problem admits at least one optimal controla*∈U.
Proof. The admissible set U is nonempty because a(t) = 0 is admissible for any u0∈[0,umax].
The positivity and boundedness analysis established earlier provides the basis for the corresponding properties of the controlled system. In particular, since u(t)≥0,
and
Furthermore,
Thus, the biological state variables remain bounded on the finite interval [0,T].
Since the admissible controls are uniformly bounded in L∞(0,T), any minimizing sequence {an} admits a subsequence converging weakly in L2(0,T) to some a*. The control bounds are preserved under weak convergence because they define a closed convex subset of L2(0,T).
Moreover, the corresponding treatment profiles
are uniformly bounded and equi-Lipschitz. By the Arzelà–Ascoli theorem, a subsequence converges uniformly to
which preserves the treatment-state constraints. Consequently, a*∈U.
The boundedness of the biological states and the regularity of the state equations imply, along a subsequence, uniform convergence of the associated state trajectories to the solution corresponding to a*.
Finally, the running cost
is continuous in the state variables and convex in a. The state-dependent terms converge under uniform convergence, while the quadratic control cost is weakly lower semicontinuous. Hence,
Therefore, a* is an optimal control.
Hamiltonian and adjoint system
Let λM, λI, λR, and λu denote the adjoint variables associated with the states M, I, R, and u, respectively.
The Hamiltonian is
On intervals where the treatment-intensity state constraint is inactive, Pontryagin’s Maximum Principle gives
Recall that
Therefore,
The adjoint system is consequently
The linear treatment-intensity cost modifies the adjoint equation associated with u through the term −C. Consequently, although the algebraic optimality condition for a remains unchanged, the resulting treatment strategy depends on C through the adjoint variable λu.
The terminal cost is
If the treatment-intensity constraint is inactive at the terminal time, the transversality conditions are
Characterization of the optimal control
On intervals where
the Hamiltonian is minimized with respect to the treatment-adjustment rate a.
The stationarity condition
gives the unconstrained candidate
Taking into account the adjustment-rate bounds, the optimal control on interior treatment-state arcs is characterized by
where Π denotes projection onto the admissible interval.
The corresponding treatment intensity satisfies
Treatment-state constraints and numerical implementation
Because u(t) is an additional state variable, the treatment-intensity bounds constitute state constraints. In particular, feasibility requires the boundary-direction conditions
On any interval where the treatment intensity remains at one of these boundaries, a(t) = 0 almost everywhere.
In the numerical implementation, we employ a viability-projected forward–backward sweep method (FBSM), using the fourth-order Runge–Kutta method for the state and adjoint equations. The treatment-adjustment rate is projected onto a feasible interval at each time step so that the resulting treatment profile respects the prescribed intensity and adjustment-rate bounds.
The treatment-intensity penalty Cu is included in both the objective evaluation and the treatment-state adjoint equation. The computed treatment strategy is assessed using convergence diagnostics, constraint violations, and comparisons of the objective functional with the no-treatment baseline.
The projection procedure enforces numerical feasibility of the treatment profile. However, it should be distinguished from a complete state-constrained Pontryagin formulation, which may require additional multipliers and boundary-arc conditions when treatment-intensity constraints are active.
5 Numerical simulations and optimal control results
The numerical simulations presented in this section are used to validate the analytical results and to provide a biological interpretation of the model dynamics. In particular, we investigate how the treatment parameter τ influences immune activity, myelin dynamics, and the transition between persistent autoimmune oscillations and immune suppression.
5.1 Time series solutions
The analytical results established that the reduced immune subsystem is governed by a critical treatment threshold τc, which separates two qualitatively distinct regimes. When τ<τc, the immune-free equilibrium E0 is unstable and the interior equilibrium Ei exists as a nonlinear center surrounded by closed periodic orbits.
These periodic orbits are a property of the reduced frozen-myelin subsystem. In the full model, where the myelin level evolves slowly, the numerical simulations exhibit recurrent oscillatory immune dynamics that are qualitatively consistent with the reduced-system behavior, but exact periodicity is not asserted.
In this regime, the full model therefore exhibits recurrent autoimmune oscillations associated with repeated inflammatory activity. In contrast, when τ>τc, the interior equilibrium Ei does not exist and the immune-free equilibrium E0 becomes globally asymptotically stable, corresponding to suppression of sustained autoimmune activity.
To illustrate these two regimes, we select representative treatment values on either side of the threshold τc and examine the corresponding time-series and phase-plane dynamics.
The baseline parameter values and associated units used throughout the simulations are listed in Table 2. Parameter values were selected from the literature whenever available and supplemented by biologically reasonable representative values when direct experimental estimates were unavailable. Consequently, the numerical simulations are intended primarily to illustrate the qualitative behavior of the model and the influence of the principal mechanisms governing immune activation, immune regulation, myelin damage, and treatment response.
The numerical results should therefore be interpreted as qualitative and mechanistic insights into the disease dynamics rather than quantitatively calibrated patient-specific predictions.
The treatment values associated with the representative myelin states are
Using the baseline parameter values, we obtain
Thus, τc is the critical treatment threshold separating persistent autoimmune activity from immune suppression, while τi and τ0 are reference treatment values corresponding to the representative myelin levels Mi and M0, respectively.
Unless otherwise stated, all numerical simulations are performed with the initial condition [(M(0),I(0),R(0),N(0)) = (19, 2, 2, 15).] This initial condition represents a state with substantial healthy myelin together with elevated autoreactive and regulatory immune cell populations, allowing the transient and asymptotic dynamics of the model to be illustrated under different treatment levels.
Representative values of τ are selected from the different regions determined by these treatment levels to illustrate the corresponding model dynamics.
The corresponding dynamics are illustrated in Fig 3, which combines the myelin, immune-cell, and phase-plane results into three panels.
Panel A shows the evolution of healthy myelin M(t) and neuronal capacity N(t). The healthy myelin level decreases over time toward the representative inflammatory myelin level, while neuronal capacity declines as damage accumulates. These dynamics illustrate the effects of persistent autoimmune activity on myelin preservation and neuronal function within the proposed model.
For the parameter values used in the simulation, the representative inflammatory state is
Here, (Ii,Ri) denotes the interior equilibrium of the reduced immune subsystem, while Mi is the corresponding representative myelin level.
Since R0>1, the immune-free equilibrium of the reduced subsystem is unstable, whereas its interior equilibrium is a nonlinear center surrounded by closed periodic orbits. This theoretical result characterizes the reduced frozen-myelin dynamics and does not, by itself, establish exact periodicity in the full model.
Panel B presents the corresponding immune-cell dynamics. The autoreactive immune-cell population I(t) exhibits recurrent peaks accompanied by oscillations in the regulatory immune-cell population R(t). This interaction produces repeated episodes of elevated and reduced immune activity, qualitatively representing recurrent inflammatory dynamics within the model.
Panel C shows the (I,R) phase-plane projection of the full-system trajectory. The recurrent oscillatory structure is qualitatively consistent with the closed-orbit behavior established for the reduced immune subsystem. However, exact periodicity is not asserted for the full system because the myelin level evolves over time.
The corresponding dynamics are illustrated in Fig 4, which presents the myelin, immune-cell, and phase-plane results in three panels.
Panel A shows the evolution of healthy myelin M(t) and neuronal capacity N(t). As autoreactive immune activity is suppressed, the healthy myelin level approaches its immune-free equilibrium value. Neuronal capacity decreases during the transient period of myelin deficiency and subsequently changes more slowly as myelin recovers.
Since neuronal recovery is not included in the model, the accumulated loss of neuronal capacity is not reversed.
Panel B presents the evolution of autoreactive immune cells I(t) and regulatory immune cells R(t). Both populations decrease toward zero, illustrating the suppression of persistent autoimmune activity under the higher treatment level.
Panel C shows the corresponding (I,R) phase-plane trajectory. The numerical trajectory approaches the immune-free state (I,R) = (0,0), consistent with the global asymptotic stability established analytically for the reduced immune subsystem when ℛ0<1.
For the (M,I,R) subsystem, the immune-free equilibrium is
The neuronal variable N(t) is not included in this equilibrium representation because it does not feed back into the myelin–immune dynamics.
The numerical simulations of the full model are qualitatively consistent with the analytical results obtained for the reduced immune subsystem. When τ<τc, recurrent oscillatory immune dynamics are observed, consistent with the nonlinear-center and closed-orbit structure of the reduced system. In contrast, when τ>τc, the numerical trajectories show decreasing autoreactive and regulatory immune-cell populations, together with recovery of the myelin level toward its immune-free equilibrium.
These findings illustrate the treatment-dependent transition between recurrent immune activity and immune suppression predicted by the reduced-system analysis, while recognizing that the global stability theorem applies specifically to the reduced immune subsystem.
Fig 5 presents the phase portraits for τ=1.335714, 1.435714, 1.485714, 1.515714, and 1.635714. The corresponding initial conditions are (I(0),R(0)) = (1.61,0.44), (1.61,0.22), (1.61,0.11), (1.61,0.165), and (2.00,0.35), respectively. For τ<τc=1.535714, the immune-free equilibrium E0 is unstable, whereas the interior equilibrium Ei(τ) is a nonlinear center surrounded by closed periodic orbits, indicating sustained oscillations of the immune cell populations.
As τ approaches the critical value τc from below, the equilibrium moves toward the R = 0 axis since Ri→0 while Ii=λ/ρ remains constant. At τ=τc, the interior equilibrium reduces to the immune-only boundary equilibrium Eb=(λ/ρ,0). For τ>τc, the interior equilibrium no longer exists, and all trajectories converge to the immune-free equilibrium E0, indicating elimination of the autoreactive immune cells.
5.2 Sensitivity analysis of the autoimmune activation number at the immune-free state
To identify the parameters that most strongly influence the stability of the immune-free state, we investigate the sensitivity of the autoimmune activation number evaluated at the immune-free myelin level. Since the immune-free state corresponds to
the associated myelin damage level is
Consequently, the autoimmune activation number at the immune-free state is given by
This threshold determines whether small perturbations in the autoreactive immune-cell population can grow near the immune-free state. Its sensitivity therefore provides insight into the parameters governing the onset of autoimmune activity.
Analytical sensitivity indices.
The normalized forward sensitivity index of R0(M0) with respect to a parameter p is defined by
For positive parameter values and R0(M0)>0, this is equivalently expressed as
This dimensionless quantity measures the relative response of the autoimmune activation number to a small relative change in the corresponding parameter.
Differentiating the expression for R0(M0) gives
Using the baseline parameter values listed in Table 2, with τ=0.5, the autoimmune activation number is R0(M0)=2.5287. The corresponding normalized sensitivity indices are summarized in Table 3.
The results in Table 3 indicate that ηmax, τ, and θ have the largest absolute sensitivity indices. In particular, the unit sensitivity of ηmax shows that the autoimmune activation number changes proportionally with the maximum immune activation rate. The negative sensitivities associated with θ and τ indicate that increasing autoreactive immune-cell removal reduces R0(M0).
In contrast, Mmax, ω, α, and δ exhibit relatively small sensitivity magnitudes at the baseline parameter set. These findings identify immune activation and autoreactive immune-cell removal as the dominant mechanisms controlling the local disease threshold.
The analytical sensitivity indices characterize the response of R0(M0) to infinitesimal parameter perturbations at the baseline values. However, they do not describe the combined effects of simultaneous parameter variations over prescribed uncertainty ranges. To complement the analytical results, we perform a global sensitivity analysis.
Global sensitivity analysis using LHS–PRCC.
Global sensitivity analysis is conducted using Latin Hypercube Sampling (LHS) combined with Partial Rank Correlation Coefficients (PRCC). This approach assesses the relative influence of model parameters on the autoimmune activation number while allowing all selected parameters to vary simultaneously.
We consider the parameter vector
where each component is independently sampled within a prescribed interval around its baseline value pi*.
Because experimentally established uncertainty ranges are unavailable for several parameters, a variation of ±10% is adopted as an exploratory uncertainty interval:
where U denotes the uniform distribution.
A total of Ns=2000 parameter combinations are generated using LHS. For each sample, the autoimmune activation number R0(M0) is calculated using its analytical expression.
The PRCC between a parameter pi and R0(M0) is computed as
where rpi and rR0 denote the residuals obtained by regressing the rank-transformed values of pi and R0(M0), respectively, on the rank-transformed values of the remaining parameters.
PRCC values range between −1 and 1. Positive coefficients indicate a positive monotonic association with R0(M0), whereas negative coefficients indicate an inverse association. The absolute magnitude reflects the strength of the partial rank correlation.
Fig 6 illustrates the relative influence of the seven parameters on the autoimmune activation number.
The maximum immune activation rate ηmax exhibits a strong positive correlation with R0(M0), whereas the removal parameters τ and θ exhibit strong negative correlations. These three parameters have the largest absolute PRCC values, indicating that they are the dominant contributors to variations in the autoimmune activation number over the investigated parameter ranges.
The parameters Mmax and ω show weaker positive and negative associations, respectively. The PRCC values corresponding to α and δ are close to zero, suggesting limited influence within the considered uncertainty intervals.
The global sensitivity results are qualitatively consistent with the analytical sensitivity indices. In particular, both approaches identify immune activation and autoreactive immune-cell removal as the principal mechanisms affecting the autoimmune activation threshold.
This agreement supports the importance of targeting immune activation and autoreactive immune-cell clearance when developing strategies aimed at reducing the autoimmune activation number.
It should be emphasized that these conclusions are conditional on the chosen uncertainty intervals and the assumption of approximately monotonic parameter-output relationships. The analysis does not establish parameter influence outside the investigated ranges.
5.3 Optimal control solution
In this section, we compute the optimal treatment strategy for the proposed MS model by numerically solving the associated optimal control problem using a forward–backward sweep algorithm. Starting from an initial admissible control, the state system is solved forward in time while the adjoint system is solved backward in time. The control is updated at each iteration using the optimality condition (34) and projected onto the admissible control set.
The iterative process continues until convergence is achieved, that is,
where ε1 is a prescribed tolerance. The resulting optimal treatment is then used to compare disease dynamics with and without intervention. The comparison focuses on the evolution of autoreactive immune cells, the treatment intensity and its rate of adjustment, as well as the instantaneous and cumulative treatment costs.
Numerical simulations show that the optimal therapy suppresses large inflammatory flares and maintains the immune activity at controlled levels while producing clinically realistic treatment profiles due to the inclusion of rate–smoothing penalties.
Fig 7 presents the optimal-control results, combining the immune-cell dynamics, cumulative objective value, treatment intensity, and treatment-adjustment rate into four panels.
Panel A compares the evolution of the autoreactive immune-cell population I(t) under untreated and optimally controlled conditions. In the absence of treatment, the system exhibits pronounced oscillations associated with immune-regulatory feedback. Under optimal control, these oscillations are substantially reduced, and the autoreactive immune-cell population generally remains at lower levels.
Temporary crossings between the controlled and uncontrolled trajectories may occur because the objective functional minimizes a weighted combination of cumulative immune burden and treatment costs rather than imposing pointwise minimization of I(t).
Panel B shows the cumulative objective value over the treatment horizon. Compared with the untreated case, the optimal strategy reduces the accumulated objective by limiting periods of elevated autoimmune activity while accounting for treatment intensity and adjustment costs. This comparison illustrates the effectiveness of the computed intervention within the proposed mathematical framework.
Panel C displays the optimal treatment intensity u(t). The computed strategy initially increases treatment intensity toward its admissible upper bound and subsequently maintains a relatively high level of intervention over much of the treatment horizon. This behavior reflects the balance between suppressing autoreactive immune activity and penalizing treatment exposure and rapid adjustments.
Panel D compares the optimal treatment-adjustment rate a(t) with the numerical derivative u˙(t) obtained from the treatment profile. Their close agreement provides a numerical consistency check for the implementation of the optimality system and the treatment dynamics u˙=a.
Taken together, the four panels demonstrate that the computed optimal-control strategy substantially reduces recurrent autoimmune activity while accounting for treatment-related costs and adjustment constraints. These findings illustrate the potential of the proposed optimal-control framework to investigate treatment-dependent immune suppression within the mathematical model.
Overall, the effectiveness of the proposed treatment strategy is assessed by comparing the controlled and uncontrolled simulations. In the absence of treatment, recurrent autoimmune activity persists, leading to sustained inflammatory episodes and progressive myelin damage. Under the optimal treatment strategy, autoreactive immune activity is substantially suppressed, the cumulative objective value is reduced, and healthy myelin is better preserved.
These numerical findings demonstrate the potential effectiveness of treatment-induced immune regulation within the proposed mathematical framework. The rate-smoothing mechanism further supports a gradual treatment profile by penalizing rapid adjustments in treatment intensity. Consequently, the optimal control framework provides a quantitative basis for evaluating treatment strategies and investigating their potential effects on autoimmune activity and disease progression.
5.3.1 Effect of optimal treatment.
To quantify the impact of the proposed optimal control strategy, we compare the controlled solution with the untreated case over the simulation horizon. Following the approach in [8], we define the cumulative autoimmune burden by
where I(t) denotes the autoreactive immune-cell population. This quantity measures the total autoimmune activity accumulated over the treatment period.
The corresponding treatment efficiency is defined by
where At is the cumulative autoimmune burden under optimal treatment and A0 is the corresponding value in the absence of treatment. Thus, larger values of F indicate a greater reduction in cumulative autoimmune activity.
Table 4 summarizes the quantitative effect of the optimal treatment. In the absence of treatment, the cumulative autoimmune burden is
whereas under optimal treatment it is reduced to
This gives
Hence, the optimal treatment reduces the cumulative autoimmune burden by approximately 97.3% relative to the untreated case.
The objective value is also substantially reduced, from
without treatment to
under optimal treatment. This reduction indicates that the controlled strategy achieves strong suppression of autoimmune activity while maintaining a relatively small treatment-adjustment cost.
The numerical trajectories in Fig 7 provide a complementary view of these results. In the untreated case, the autoreactive immune-cell population exhibits recurrent large-amplitude inflammatory peaks. Under optimal treatment, these oscillations are rapidly suppressed and the immune-cell population remains close to zero for most of the simulation interval.
The optimal treatment profile initially increases from zero, reaches its maximum level during the early phase of the intervention, and then gradually decreases as autoimmune activity is brought under control. The corresponding treatment-adjustment rate a(t) is initially positive, allowing treatment intensity to increase, and later becomes negative once treatment intensity begins to decrease. This behavior is consistent with the relation
The close agreement between the optimal rate a(t) obtained from Pontryagin’s Maximum Principle and the numerical derivative du/dt provides an additional consistency check for the numerical solution.
Since autoreactive immune cells are the principal drivers of myelin destruction in the model, the substantial reduction in cumulative autoimmune activity is expected to reduce immune-mediated myelin damage. Overall, the results demonstrate that the proposed optimal-control strategy can strongly suppress autoimmune activity while avoiding unnecessarily large or abrupt changes in treatment intensity.
6 Conclusion
In this work, we developed and analyzed a mathematical model describing the interaction between autoreactive immune cells, regulatory immune responses, and myelin dynamics in multiple sclerosis. The model incorporates damage-dependent immune activation together with treatment-induced modulation, providing a biologically consistent and analytically tractable framework.
We established well-posedness of the system and derived a reduced immune subsystem by exploiting the separation of time scales between myelin and immune dynamics. This reduction enabled a rigorous analysis of equilibrium structure and stability properties. In particular, we identified the immune-free, boundary, and interior equilibria and characterized their existence and stability in terms of the system parameters.
A key contribution of this work is the identification of a treatment-dependent threshold structure governing the qualitative behavior of the reduced immune subsystem. While the autoimmune activation number R0 determines the local stability of the immune-free equilibrium, the analysis also reveals distinct dynamical regimes on either side of the critical threshold R0=1.
In particular, when R0>1, the interior equilibrium exists as a nonlinear center and is surrounded by closed periodic orbits, corresponding to sustained autoimmune oscillations within the reduced immune subsystem. When R0<1, the interior equilibrium no longer exists and the immune-free equilibrium is globally asymptotically stable in the biologically relevant domain of the reduced system.
At the critical threshold R0=1, the reduced system possesses a continuum of boundary equilibria, with Eb=(λ/ρ,0) representing the limiting position of the interior equilibrium. These findings identify a critical treatment-dependent threshold separating persistent autoimmune activity from immune suppression and demonstrate how treatment alters the qualitative structure of the reduced immune dynamics.
Numerical simulations of the full model are qualitatively consistent with the analytical findings for the reduced immune subsystem and reveal distinct dynamical regimes. For R0>1, the system exhibits recurrent inflammatory oscillations associated with progressive myelin damage. For sufficiently strong treatment, immune activity is suppressed and the system approaches the immune-free state.
Both analytical and global sensitivity analyses identify the maximum immune activation rate ηmax, the treatment-induced removal rate τ, and the natural removal rate θ as the most influential parameters affecting the autoimmune activation number. The agreement between normalized sensitivity indices and LHS–PRCC results highlights the importance of immune activation and autoreactive immune-cell removal in governing the disease threshold. These findings are conditional on the parameter ranges considered in the global sensitivity analysis.
An optimal control framework was formulated to identify treatment strategies that minimize disease burden while maintaining realistic therapeutic profiles. The numerical results indicate that the computed treatment strategies can reduce cumulative autoimmune activity while accounting for treatment intensity and adjustment costs in the objective functional.
Overall, the proposed framework integrates mechanistic modeling, threshold and stability analysis, local and global sensitivity analysis, and optimal control. The results provide insight into the role of immune activation, regulatory responses, and treatment in determining autoimmune persistence and suppression, while highlighting the importance of treatment-dependent changes in the qualitative structure of the reduced immune subsystem.
Parameter estimation, model calibration, and validation were not pursued in the present study because suitable clinical or experimental datasets containing all model variables are not currently available. The primary objective was to develop and analyze a mechanistic mathematical framework capable of capturing qualitative features of immune dynamics in multiple sclerosis. Consequently, the numerical findings should be interpreted as mechanistic predictions rather than patient-specific clinical outcomes.
Moreover, the present study is based on a deterministic ordinary differential equation framework and therefore does not explicitly account for random fluctuations associated with immune heterogeneity, environmental variability, or uncertainty in treatment response.
In addition, the neuronal equation represents cumulative neuronal loss and does not include neuronal repair or feedback from neuronal damage to the immune dynamics. Incorporating neuronal recovery and bidirectional neuroimmune feedback would therefore be a relevant extension of the present framework.
Future work will focus on parameter estimation and validation using suitable clinical datasets, together with extensions incorporating additional immune pathways, patient heterogeneity, and time delays associated with immune activation, regulatory responses, demyelination, remyelination, and neuronal damage. Stochastic formulations may also be investigated to examine noise-induced transitions, relapse–remission patterns, and uncertainty in treatment outcomes. Such developments may improve the predictive capability and biological relevance of the proposed framework and reveal additional dynamical behaviors, including delay-induced oscillations and bifurcations.
These directions are motivated by recent advances in stochastic disease modeling, delay differential equations, Hopf bifurcation analysis, and optimal control of demyelinating diseases [9–12].