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: a.
The model represents an internal autoimmune process with no external infection, environmental transmission, or migration. b. Autoreactive immune cells become activated in response to myelin damage in the central nervous system. c. Regulatory immune cells suppress the proliferation and inflammatory activity of autoreactive immune cells. d.
Myelin is subject to intrinsic turnover and immune-mediated degradation. e. Remyelination and demyelination occur on a slower physiological timescale than immune cell activation and suppression. f. Homogeneous mixing is assumed, meaning that cellular interactions occur uniformly within the system.
The homogeneous mixing assumption is a standard approximation in compartmental models and implies that interactions among immune cells and myelin are represented by average population densities rather than explicit spatial distributions. Although multiple sclerosis exhibits spatial heterogeneity and localized lesions, the present study focuses on the overall immune-myelin dynamics and treatment effects at the population level. g. Myelin has a theoretical upper bound M max corresponding to a fully healthy state, which is not exceeded in biological conditions.
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. T reg cells) R ( t ). Autoreactive immune activation is assumed to increase with the level of myelin damage.
We define the damage variable by D = M max − M , where M max 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 η ( D ) = η max D ω + D , η max > 0 , ω > 0 , 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 β I M . Autoreactive immune cells I ( t ) proliferate at rate η ( D ) I , are suppressed by regulatory cells at rate μ R I , and are removed by natural death and treatment a t rate ( θ + τ ) I , where θ denotes natural cell death and τ reflects therapeutic intervention. Regulatory immune cells R ( t ) expand proportionally to autoreactive cell presence at rate ρ I R and decay naturally at rate λ R .
The functional capacity of neurons N ( t ) deteriorates when myelin falls below its maximal healthy level M max . 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 . 10.1371/journal.pone.0349131.t001 Table 1 Description of model parameters.
Parameter Description α Natural regeneration rate of healthy myelin β Rate of immune-mediated myelin degradation δ Natural degradation rate of myelin η max Maximum activation rate of autoreactive immune cells ω Half-saturation constant in the activation function η ( D ) μ Suppression rate of autoreactive immune cells by regulatory cells θ Natural death rate of autoreactive immune cells τ Treatment-induced removal rate of autoreactive immune cells ρ Proliferation rate of regulatory immune cells in the presence of inflammation λ Natural decay rate of regulatory immune cells M max Maximum attainable level of healthy myelin ε Small parameter representing the slower myelin dynamics relative to immune responses γ Rate of neuronal loss due to myelin deficiency 10.1371/journal.pone.0349131.g001 Fig 1 Flowchart of the proposed multiple sclerosis model. The damage variable is defined in the text by D = M max − M . Healthy myelin M is regenerated at rate α , degraded naturally at rate δ M , and damaged by autoreactive immune cells at rate β I M .
Myelin damage activates autoreactive immune cells through the damage-dependent term η ( D ) I . Regulatory immune cells suppress autoreactive immune cells at rate μ R I , expand at rate ρ I R , and decay at rate λ R . Neuronal capacity N declines due to myelin damage at rate γ D N .
The model resulting from these assumptions is given by the following system of ordinary differential equations: d M d t = ε ( α − β I M − δ M ) , (1) d I d t = ( η ( D ) − μ R − ( θ + τ ) ) I , (2) d R d t = ( ρ I − λ ) R , (3) d N d t = − ε γ D N (4) 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. 10.1371/journal.pone.0349131.g002 Fig 2 Schematic illustration of the treatment-dependent stability and threshold structure of the reduced immune subsystem.
The critical treatment threshold τ c separates two qualitatively distinct dynamical regimes. For τ < τ c (equivalently, ℛ 0 > 1 ), the interior equilibrium E i exists as a nonlinear center and is surrounded by closed periodic orbits, representing sustained autoimmune oscillations. At the critical threshold τ = τ c (equivalently, ℛ 0 = 1 ), the system possesses the continuum of boundary equilibria ℰ c = { ( I , 0 ) : I ≥ 0 } , and E b = ( λ / ρ , 0 ) is the limiting position of E i as τ → τ c − .
For τ > τ c (equivalently, ℛ 0 < 1 ), the interior equilibrium no longer exists, and the immune-free equilibrium E 0 is globally asymptotically stable, with all trajectories converging to E 0 . 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 𝐅 ( t , 𝐗 ) , 𝐗 ( t ) = ( M ( t ) , I ( t ) , R ( t ) , N ( t ) ) ⊤ . 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 = M max − 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 ≤ M max . Therefore, by the Picard–Lindelöf (Cauchy–Lipschitz) theorem, for any initial condition 𝐗 ( 0 ) ∈ ℝ + 4 , there exists a unique local solution 𝐗 ( t ) = ( M ( t ) , I ( t ) , R ( t ) , N ( t ) ) ⊤ defined on some interval [0, T 1 ), where T 1 > 0 is the maximal time of existence.
2.4 Nonnegativity of Solutions Theorem 1. Assume the initial data are nonnegative, M ( 0 ) ≥ 0 , I ( 0 ) ≥ 0 , R ( 0 ) ≥ 0 , N ( 0 ) ≥ 0 . Then the corresponding solution satisfies M ( t ) ≥ 0 , I ( t ) ≥ 0 , R ( t ) ≥ 0 , N ( t ) ≥ 0 for all t ∈ [ 0 , T ) .
Proof. We show that the nonnegative orthant is forward invariant. On the boundary M = 0, M ˙ | M = 0 = ε α > 0 .
Thus, M ( t ) cannot cross into negative values. The I -equation has the form I ˙ = ( η ( D ) − μ R − ( θ + τ ) ) I , so I ( t ) = I ( 0 ) exp ( ∫ 0 t ( η ( D ( s ) ) − μ R ( s ) − ( θ + τ ) ) d s ) ≥ 0 . Similarly, R ( t ) = R ( 0 ) exp ( ∫ 0 t ( ρ I ( s ) − λ ) d s ) ≥ 0 .
Finally, N ˙ = − ε γ D N , so N ( t ) = N ( 0 ) exp ( − ε γ ∫ 0 t D ( s ) d s ) ≥ 0 . 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, M ˙ = ε ( α − β I M − δ M ) ≤ ε ( α − δ M ) , since I ≥ 0 . By the comparison principle, M ( t ) ≤ M ( 0 ) e − ε δ t + α δ ( 1 − e − ε δ t ) . Hence 0 ≤ M ( t ) ≤ max { M ( 0 ) , α δ } .
Next, observe that η ( D ) = η max D ω + D , D ≥ 0 , satisfies 0 ≤ η ( D ) ≤ η max . Therefore, I ˙ = ( η ( D ) − ( θ + τ ) − μ R ) I ≤ ( η max − ( θ + τ ) − μ R ) I , while R ˙ = ( ρ I − λ ) R . We distinguish two cases.
Case 1: η max ≤ θ + τ . In this case, I ˙ ≤ − ( θ + τ − η max ) I , and therefore Let k = θ + τ − η max > 0 . Then I ( t ) ≤ I ( 0 ) e − k t .
Hence I ( t ) is bounded and converges exponentially to zero. For R ( t ), we write R ˙ = ( ρ I ( t ) − λ ) R . Solving this scalar equation gives R ( t ) = R ( 0 ) exp ( ∫ 0 t ( ρ I ( s ) − λ ) d s ) .
Using the bound on I ( t ), we obtain ∫ 0 t I ( s ) d s ≤ I ( 0 ) ∫ 0 t e − k s d s = I ( 0 ) k ( 1 − e − k t ) ≤ I ( 0 ) k . Therefore, R ( t ) ≤ R ( 0 ) exp ( ρ I ( 0 ) k ) . Thus R ( t ) remains uniformly bounded on [ 0 , ∞ ) .
Consequently, both immune populations are bounded. Case 2: η max > θ + τ . Consider the Lotka–Volterra comparison system I ˙ = ( η max − ( θ + τ ) − μ R ) I , R ˙ = ( ρ I − λ ) R .
This system admits the positive equilibrium I i = λ ρ , R i = η max − ( θ + τ ) μ . Moreover, V ( I , R ) = ρ I − λ ln I + μ R − ( η max − ( θ + τ ) ) ln R is a first integral satisfying d V d t = 0 . Furthermore, V ( I , R ) → ∞ as I → 0 + , R → 0 + , I → ∞ , R → ∞ .
Hence every level set of V is compact in the positive quadrant. Let V 0 = V ( I ( 0 ) , R ( 0 ) ) and define Γ V 0 = { ( I , R ) ∈ ℝ + 2 : V ( I , R ) ≤ V 0 } . Since Γ V 0 is compact and the function Φ ( I , R ) = I + R is continuous, there exists Φ m a x = max ( I , R ) ∈ Γ V 0 ( I + R ) .
Therefore, 0 ≤ I ( t ) + R ( t ) ≤ Φ m a x , t ≥ 0 , where Φ m a x = Φ m a x ( I ( 0 ) , R ( 0 ) , ρ , λ , μ , η max , θ , τ ) > 0 . Since the original immune subsystem is dominated by the above comparison system, its trajectories remain bounded as well. Finally, for the neuronal capacity variable, N ˙ = − ε γ D N , with D ≥ 0 .
Therefore N ˙ ≤ 0 , which implies 0 ≤ N ( t ) ≤ N ( 0 ) , t ≥ 0 . 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 Ω = { ( M , I , R , N ) ∈ ℝ + 4 : 0 ≤ M ≤ M max * , 0 ≤ I + R ≤ Φ * , 0 ≤ N ≤ N ( 0 ) } , where M max * = max { M ( 0 ) , α δ } , and Φ * = max { Φ max , I ( 0 ) + R ( 0 ) exp ( ρ I ( 0 ) θ + τ − η max ) } .
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, 0 ≤ M ( t ) ≤ M max * , 0 ≤ I ( t ) ≤ I ( t ) + R ( t ) ≤ Φ * .
Hence, from M ˙ = ε ( α − β I M − δ M ) , we obtain | M ( t ) − M ( 0 ) | ≤ ε ∫ 0 t ( α + β I ( s ) M ( s ) + δ M ( s ) ) d s ≤ ε t [ α + ( β Φ * + δ ) M max * ] ≤ ε T [ α + ( β Φ * + δ ) M max * ] . (5) Therefore, M ( t ) = M ( 0 ) + O ( ε ) uniformly for 0 ≤ t ≤ T . Since D ( t ) = M max − M ( t ) , it follows that D ( t ) = D ( 0 ) + O ( ε ) .
Moreover, η ( D ) = η max D ω + D , η ′ ( D ) = η max ω ( ω + D ) 2 ≤ η max ω , D ≥ 0 . Thus, using (5) , | η ( D ( t ) ) − η ( D ( 0 ) ) | ≤ η max ω | D ( t ) − D ( 0 ) | = η max ω | M ( t ) − M ( 0 ) | ≤ ε η max ω T [ α + ( β Φ * + δ ) M max * ] . (6) To compare the full and frozen-myelin immune dynamics, define X ε ( t ) = ( I ε ( t ) R ε ( t ) ) , X 0 ( t ) = ( I 0 ( t ) R 0 ( t ) ) , and F ( X , M ) = ( [ η ( M max − M ) − μ R − θ − τ ] I ( ρ I − λ ) R ) .
Then X ˙ ε = F ( X ε , M ( t ) ) , X ˙ 0 = F ( X 0 , M ( 0 ) ) , with X ε ( 0 ) = X 0 ( 0 ) . For the same immune state X = ( I , R ) ⊤ , F ( X , M ( t ) ) − F ( X , M ( 0 ) ) = ( [ η ( D ( t ) ) − η ( D ( 0 ) ) ] I 0 ) . Using I ≤ Φ * and (6) , ‖ F ( X , M ( t ) ) − F ( X , M ( 0 ) ) ‖ ≤ Φ * | η ( D ( t ) ) − η ( D ( 0 ) ) | ≤ ε K T , (7) where K T = Φ * η max ω T [ α + ( β Φ * + δ ) M max * ] .
Since F is locally Lipschitz in X on the bounded invariant region Ω , there exists L T > 0 such that ‖ X ε ( t ) − X 0 ( t ) ‖ ≤ L T ∫ 0 t ‖ X ε ( s ) − X 0 ( s ) ‖ d s + ε K T t . Grönwall’s inequality therefore gives ‖ X ε ( t ) − X 0 ( t ) ‖ ≤ ε K T L T ( e L T t − 1 ) , and consequently sup 0 ≤ t ≤ T ‖ X ε ( t ) − X 0 ( t ) ‖ = O ( ε ) . 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 D 0 = M max − M ( 0 ) . The corresponding fast layer problem is d I d t = ( η ( D 0 ) − μ R − θ − τ ) I , d R d t = ( ρ I − λ ) R . (8) For notational simplicity, the prescribed damage level D 0 will henceforth be denoted by D .
Therefore, (8) becomes d I d t = ( η ( D ) − μ R − θ − τ ) I , d R d t = ( ρ I − λ ) R , 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 η ( D ) > θ + τ , 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 I ( η ( D ) − μ R − θ − τ ) = 0 , and R ( ρ I − λ ) = 0 .
The immune-free equilibrium is E 0 = ( 0 , 0 ) . For a positive interior equilibrium E i , I i = λ ρ , and R i = η ( D ) − θ − τ μ . Hence, the immune subsystem admits the interior equilibrium E i = ( λ ρ , η ( D ) − θ − τ μ ) , provided that η ( D ) > θ + τ .
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 η ( D ) = θ + τ , the interior equilibrium reduces to the immune-only boundary equilibrium E b = ( λ ρ , 0 ) . 3.2 Local stability analysis Let F 1 ( I , R ) = ( η ( D ) − μ R − θ − τ ) I , (9) and F 2 ( I , R ) = ( ρ I − λ ) R .
(10) The Jacobian matrix of system (8) is J ( I , R ) = ( η ( D ) − μ R − θ − τ − μ I ρ R ρ I − λ ) . (11) 3.2.1 Immune-free equilibrium. The immune-free equilibrium of the immune subsystem is E 0 = ( 0 , 0 ) .
The Jacobian matrix evaluated at E 0 is J ( E 0 ) = ( η ( D ) − θ − τ 0 0 − λ ) . (12) Hence, the corresponding eigenvalues are ξ 1 = η ( D ) − θ − τ , ξ 2 = − λ . (13) Therefore, the immune-free equilibrium is locally asymptotically stable whenever η ( D ) < θ + τ , (14) and unstable whenever η ( D ) > θ + τ .
(15) Motivated by this threshold condition, we define the damage-dependent autoimmune activation number R 0 = η ( D ) θ + τ = η max θ + τ D ω + D . (16) Using D = M max − M , this may also be written as R 0 = η max θ + τ M max − M ω + M max − M . Consequently, R 0 < 1 ⟺ η ( D ) < θ + τ , whereas R 0 > 1 ⟺ η ( D ) > θ + τ .
The quantity R 0 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 R 0 < 1 , immune removal dominates damage-induced activation and autoreactive immune activity declines.
In contrast, when R 0 > 1 , damage-induced activation exceeds immune removal, causing autoreactive immune activity to increase. Consequently, R 0 = 1 represents the threshold separating immune control from immune activation. The threshold condition may also be expressed in terms of the treatment parameter.
Defining τ c = η ( D ) − θ , (17) we obtain τ > τ c ⟺ R 0 < 1 , and τ < τ c ⟺ R 0 > 1 . 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 η ( D ) > θ + τ , (18) so that the interior equilibrium E i = ( λ ρ , η ( D ) − θ − τ μ ) (19) exists.
Evaluating the Jacobian at E i gives J ( E i ) = ( 0 − μ λ / ρ ρ ( η ( D ) − θ − τ ) / μ 0 ) . (20) The characteristic equation is ξ 2 + λ ( η ( D ) − θ − τ ) = 0 . (21) Thus, ξ 1 , 2 = ± i λ ( η ( D ) − θ − τ ) .
(22) 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 η ( D ) > θ + τ .
Define a = η ( D ) − θ − τ > 0 . Then the immune subsystem becomes d I d t = ( a − μ R ) I , d R d t = ( ρ I − λ ) R . (23) The interior equilibrium is E i = ( λ ρ , a μ ) .
We now show that E i is a nonlinear center. Consider the function H ( I , R ) = ρ I − λ ln I + μ R − a ln R , I > 0 , R > 0 . (24) Differentiating H along solutions gives d H d t = ( ρ − λ I ) d I d t + ( μ − a R ) d R d t .
Substituting the system equations yields d H d t = ( ρ − λ I ) ( a − μ R ) I + ( μ − a R ) ( ρ I − λ ) R . Since ( ρ − λ I ) I = ρ I − λ , and ( μ − a R ) R = μ R − a , we obtain d H d t = ( ρ I − λ ) ( a − μ R ) + ( μ R − a ) ( ρ I − λ ) . Using μ R − a = − ( a − μ R ) , it follows that d H d t = ( ρ I − λ ) ( a − μ R ) − ( a − μ R ) ( ρ I − λ ) = 0 .
Hence, H ( I , R ) = constant along every solution trajectory. Next, compute the first derivatives: H I = ρ − λ I , H R = μ − a R . Therefore, H I = 0 ⟺ I = λ ρ , and H R = 0 ⟺ R = a μ .
Thus, the only critical point of H in the positive quadrant is E i . The second derivatives are H I I = λ I 2 > 0 , H R R = a R 2 > 0 , H I R = 0 . Hence, the Hessian matrix is ∇ 2 H = ( λ / I 2 0 0 a / R 2 ) , which is positive definite throughout the positive quadrant.
Therefore, H has a strict local minimum at E i . Furthermore, H ( I , R ) → ∞ as I → 0 + , R → 0 + , I → ∞ , R → ∞ . Consequently, every level set H ( I , R ) = C , C > H ( E i ) , sufficiently close to E i is a simple closed curve surrounding E i .
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 E i are closed orbits surrounding E i . We conclude that the interior equilibrium E i = ( λ ρ , η ( D ) − θ − τ μ ) is a nonlinear center.
3.2.3 Immune-only boundary equilibrium. At the critical value η ( D ) = θ + τ , the interior equilibrium reduces to the immune-only boundary equilibrium E b = ( λ ρ , 0 ) . Evaluating the Jacobian at E b yields J ( E b ) = ( 0 − μ λ / ρ 0 0 ) , whose characteristic equation is ξ 2 = 0 .
Hence, the point E b = ( λ ρ , 0 ) possesses a double zero eigenvalue, and linearization is therefore inconclusive. This is consistent with the degenerate critical structure described above at η ( D ) = θ + τ . 3.3 Global stability of the immune-free equilibrium Theorem 3.
If R 0 < 1 , then the immune-free equilibrium E 0 = ( 0 , 0 ) of system (8) is globally asymptotically stable in ℝ + 2 . Proof. The condition R 0 < 1 is equivalent to η ( D ) < θ + τ .
Let c = θ + τ − η ( D ) > 0 . Consider the Lyapunov function V ( I , R ) = I + μ ρ R . Clearly, V ( I , R ) ≥ 0 for all ( I , R ) ∈ ℝ + 2 , and V ( I , R ) = 0 ⟺ ( I , R ) = ( 0 , 0 ) .
Differentiating V along solutions gives d V d t = d I d t + μ ρ d R d t . Using system (8) , we obtain d V d t = ( η ( D ) − μ R − θ − τ ) I + μ ρ ( ρ I − λ ) R . Expanding, d V d t = ( η ( D ) − θ − τ ) I − μ I R + μ I R − μ λ ρ R .
Thus, d V d t = ( η ( D ) − θ − τ ) I − μ λ ρ R . Since η ( D ) − θ − τ = − c < 0 , we get d V d t = − c I − μ λ ρ R ≤ 0 . Moreover.
Therefore, by Lyapunov’s direct method, the immune-free equilibrium E 0 = (0,0) is globally asymptotically stable in ℝ + 2 . 3.4 Treatment-induced stability and long-term dynamics The threshold quantity R 0 ( M ) = η ( M max − M ) θ + τ depends on the current myelin level M . Since D = M max − M , and η ( D ) is an increasing function of the damage level D , it follows that R 0 ( M ) is a continuous decreasing function of M .
Consequently, there exists a unique critical myelin level M c satisfying R 0 ( M c ) = 1 . 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 E 0 is globally asymptotically stable and the immune populations converge to zero.
In this case, the slow myelin equation reduces to M ˙ = α − δ M , yielding the disease-free myelin level M 0 = α δ . When the treatment level is below τ c , the interior equilibrium E i exists and the immune subsystem exhibits sustained oscillatory dynamics around the representative autoreactive burden I i = λ ρ . Substituting this value into the slow myelin equation gives M ˙ = α − δ M − β M λ ρ , which yields the inflammatory myelin level M i = α δ + β λ / ρ .
Since δ + β λ ρ > δ , we obtain M i < M 0 . For the parameter regime considered in this work, the disease-free state satisfies R 0 ( M 0 ) < 1 , whereas the inflammatory state satisfies R 0 ( M i ) > 1 . Since R 0 ( M ) is continuous and strictly decreasing in M , the unique threshold level M c must lie between M i and M 0 .
Therefore, 0 ≤ M i < M c < M 0 ≤ M max . Defining the corresponding damage levels by D j = M max − M j , j = i , c , 0 , gives D i > D c > D 0 . Associated with these damage levels are the treatment values τ i = η ( D i ) − θ , τ c = η ( D c ) − θ , τ 0 = η ( D 0 ) − θ .
Since η ( D ) is increasing in D , it follows that τ i > τ c > τ 0 . 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 E 0 is unstable, while the interior equilibrium E i 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 ℰ c = { ( I , 0 ) : I ≥ 0 } , with E b = ( λ ρ , 0 ) being the limiting position of E i as τ → τ c − .
For treatment levels above τ c , the interior equilibrium no longer exists and the immune-free equilibrium E 0 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 ) , represen t ing the rate of treatment adjustment. Accordingly, u ( t ) is treated as an additional state variable, while a ( t ) is the con t rol to be optimized.
Positive values of a ( t ) correspond to treatment escalation, whereas negative values represent treatment reduction. Since a ∈ L ∞ ( 0 , T ) and u ( t ) = u 0 + ∫ 0 t a ( s ) d s , the treatment intensity satisfies u ∈ W 1 , ∞ ( 0 , T ) and is therefore Lipschitz continuous. Controlled system and objective functional The controlled system is given by M ˙ ( t ) = ε ( α − β I ( t ) M ( t ) − δ M ( t ) ) , I ˙ ( t ) = ( η ( D ( M ( t ) ) ) − μ R ( t ) − θ − u ( t ) ) I ( t ) , R ˙ ( t ) = ( ρ I ( t ) − λ ) R ( t ) , u ˙ ( t ) = a ( t ) , (25) subject to M ( 0 ) = M 0 , I ( 0 ) = I 0 , R ( 0 ) = R 0 , u ( 0 ) = u 0 .
(26) The objective functional is defined by J ( a ) = ∫ 0 T [ A I ( t ) + B a ( t ) 2 + C u ( t ) ] d t + q I I ( T ) , (27) where A > 0 , B > 0 , C > 0 , q I ≥ 0 . The term AI ( t ) penalizes autoreactive immune-cell burden, while Ba ( t ) 2 penalizes rapid changes in treatment intensity. The terminal term q I I ( 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 U = { a ∈ L ∞ ( 0 , T ) : a min ≤ a ( t ) ≤ a max a.e. , 0 ≤ u 0 + ∫ 0 t a ( s ) d s ≤ u max , t ∈ [ 0 , T ] } , (28) where a min ≤ 0 ≤ a max , 0 ≤ u 0 ≤ u max .
The optimal control problem consists of finding a * ∈ U satisfying J ( a * ) = min a ∈ U J ( a ) . (29) Existence of an optimal control Theorem 4. Assume that the biological parameters are positive, the initial conditions are nonnegative, and 0 ≤ M 0 ≤ M max , α δ ≤ M max .
Then the optimal control problem admits at least one optimal control a * ∈ U . Proof. The admissible set U is nonempty because a ( t ) = 0 is admissible for any u 0 ∈ [ 0 , u max ] .
The positivity and boundedness analysis established earlier provides the basis for the corresponding properties of the controlled system. In particular, since u ( t ) ≥ 0 , 0 ≤ M ( t ) ≤ max { M 0 , α δ } ≤ M max , and 0 ≤ I ( t ) ≤ I 0 e η max T . Furthermore, 0 ≤ R ( t ) ≤ R 0 exp ( ρ I 0 T e η max T ) .
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 { a n } admits a subsequence converging weakly in L 2 (0, T ) to some a * . The control bounds are preserved under weak convergence because they define a closed convex subset of L 2 (0, T ).
Moreover, the corresponding treatment profiles u n ( t ) = u 0 + ∫ 0 t a n ( s ) d s are uniformly bounded and equi-Lipschitz. By the Arzelà–Ascoli theorem, a subsequence converges uniformly to u * ( t ) = u 0 + ∫ 0 t a * ( s ) d s , 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 L ( I , u , a ) = A I + B a 2 + C u 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, J ( a * ) ≤ lim inf n → ∞ J ( a n ) = inf a ∈ U J ( a ) . 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 H = A I + B a 2 + C u + λ M ε ( α − β I M − δ M ) + λ I ( η ( D ( M ) ) − μ R − θ − u ) I + λ R ( ρ I − λ ) R + λ u a . (30) On intervals where the treatment-intensity state constraint is inactive, Pontryagin’s Maximum Principle gives λ ˙ j = − ∂ H ∂ x j , x j ∈ { M , I , R , u } . Recall that D ( M ) = M max − M , η ( D ) = η max D ω + D .
Therefore, η M ( M ) : = d η ( D ( M ) ) d M = − η max ω ( ω + M max − M ) 2 . (31) The adjoint system is consequently λ ˙ M = ε λ M ( β I + δ ) − λ I I η M ( M ) , λ ˙ I = − A + ε β M λ M − λ I ( η ( D ( M ) ) − μ R − θ − u ) − ρ R λ R , λ ˙ R = μ I λ I − ( ρ I − λ ) λ R , λ ˙ u = λ I I − C . (32) 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 Φ = q I I ( T ) . If the treatment-intensity constraint is inactive at the terminal time, the transversality conditions are λ M ( T ) = 0 , λ I ( T ) = q I , λ R ( T ) = 0 , λ u ( T ) = 0 .
(33) Characterization of the optimal control On intervals where 0 < u * ( t ) < u max , the Hamiltonian is minimized with respect to the treatment-adjustment rate a . The stationarity condition ∂ H ∂ a = 2 B a + λ u = 0 gives the unconstrained candidate a ^ ( t ) = − λ u ( t ) 2 B . Taking into account the adjustment-rate bounds, the optimal control on interior treatment-state arcs is characterized by a * ( t ) = Π [ a min , a max ] ( − λ u ( t ) 2 B ) , (34) where Π denotes projection onto the admissible interval.
The corresponding treatment intensity satisfies u * ( t ) = u 0 + ∫ 0 t a * ( s ) d s . (35) 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 u ( t ) = 0 ⟹ a ( t ) ≥ 0 , u ( t ) = u max ⟹ a ( t ) ≤ 0 .
(36) 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 E 0 is unstable and the interior equilibrium E i 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 E i does not exist and the immune-free equilibrium E 0 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. 10.1371/journal.pone.0349131.t002 Table 2 Baseline parameter values and units used in the numerical simulations and optimal control experiments.
Parameter Value Units ε 0.05 dimensionless α 1.0 myelin·time –1 β 0.1 cell − 1 · time –1 δ 0.3 time –1 γ 10 –5 time –1 μ 0.5 cell − 1 · time –1 θ 0.25 time –1 τ variable time –1 ρ 0.25 cell − 1 · time –1 λ 0.35 time –1 M max 40 myelin η max 2.0 time –1 ω 2 myelin A 100 cost · cell − 1 · time − 1 B 3 cost · time 3 q I 500 cost · cell − 1 u min 0.0 time –1 u max 1.0 time –1 a min − 1 time –2 a max 1 time –2 The treatment values associated with the representative myelin states are τ c = η ( D c ) − θ , τ i = η ( D i ) − θ , τ 0 = η ( D 0 ) − θ . Using the baseline parameter values, we obtain τ i = 1.547235 , τ i > τ c > τ 0 , τ 0 = 1.535714 . 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 M i and M 0 , 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. Case 1: τ = 0 ( τ < τ c ) .
This case corresponds to the autoimmune regime in which R 0 > 1 . The corresponding dynamics are illustrated in Fig 3 , which combines the myelin, immune-cell, and phase-plane results into three panels. 10.1371/journal.pone.0349131.g003 Fig 3 Immune and myelin dynamics for ℛ 0 > 1 .
Time-series dynamics for τ = 0 , corresponding to ℛ 0 = 7.1429 > 1 . (A) Evolution of healthy myelin M ( t ) and neuronal capacity N ( t ). (B) Evolution of autoreactive immune cells I ( t ) and regulatory immune cells R ( t ).
(C) Phase-plane trajectory of autoreac t ive immune cells I ( t ) and regula t ory immune cells R ( t ). 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 ( M i , I i , R i ) = ( 2.272727 , 1.400000 , 3.094470 ) . Here, ( I i , R i ) denotes the interior equilibrium of the reduced immune subsystem, while M i is the corresponding representative myelin level.
Since R 0 > 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. Case 2: τ = 2 ( τ > τ c ) .
This case corresponds to the immune-suppression regime in which ℛ 0 < 1 . The corresponding dynamics are illustrated in Fig 4 , which presents the myelin, immune-cell, and phase-plane results in three panels. 10.1371/journal.pone.0349131.g004 Fig 4 Immune and myelin dynamics for ℛ 0 < 1 .
Time-series dynamics for τ = 2 , corresponding to ℛ 0 = 0.7937 < 1 . (A) Evolution of healthy myelin M ( t ) and neuronal capacity N ( t ). (B) Evolution of autoreactive immune cells I ( t ) and regulatory immune cells R ( t ).
(C) Phase-plane trajectory of autoreac t ive immune cells I ( t ) and regula t ory immune cells R ( t ). 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 ( M 0 , I 0 , R 0 ) = ( α δ , 0 , 0 ) = ( 3.333333 , 0 , 0 ) . 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 E 0 is unstable, whereas the interior equilibrium E i ( τ ) 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 R i → 0 while I i = λ / ρ remains constant. At τ = τ c , the interior equilibrium reduces to the immune-only boundary equilibrium E b = ( λ / ρ , 0 ) . For τ > τ c , the interior equilibrium no longer exists, and all trajectories converge to the immune-free equilibrium E 0 , indicating elimination of the autoreactive immune cells.
10.1371/journal.pone.0349131.g005 Fig 5 Phase portraits for different treatment levels τ . For τ < τ c , the full-model trajectories exhibit recurrent oscillatory behavior consistent with the nonlinear-center and closed-orbit structure of the reduced immune subsystem. As τ → τ c − , the interior equilibrium satisfies R i → 0 and approaches E b = ( λ / ρ , 0 ) , which lies on the critical equilibrium set ℰ c = { ( I , 0 ) : I ≥ 0 } .
For τ > τ c , the immune variables converge toward the immune-free equilibrium E 0 , indicating suppression of sustained autoreactive immune activity. 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 M 0 = α δ , the associated myelin damage level is D 0 = M max − M 0 = M max − α δ .
Consequently, the autoimmune activation number at the immune-free state is given by R 0 ( M 0 ) = η ( D 0 ) θ + τ = η max θ + τ M max − α δ ω + M max − α δ . 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 R 0 ( M 0 ) with respect to a parameter p is defined by Υ p R 0 = p R 0 ( M 0 ) ∂ R 0 ( M 0 ) ∂ p . For positive parameter values and R 0 ( M 0 ) > 0 , this is equivalently expressed as Υ p R 0 = ∂ ln R 0 ( M 0 ) ∂ ln p .
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 R 0 ( M 0 ) gives Υ η max R 0 = 1 , Υ θ R 0 = − θ θ + τ , Υ τ R 0 = − τ θ + τ , Υ M max R 0 = M max δ 2 ω ( δ ( M max + ω ) − α ) ( M max δ − α ) , Υ α R 0 = − α δ ω ( δ ( M max + ω ) − α ) ( M max δ − α ) , Υ δ R 0 = α δ ω ( δ ( M max + ω ) − α ) ( M max δ − α ) , Υ ω R 0 = − δ ω δ ( M max + ω ) − α . (37) Using the baseline parameter values listed in Table 2 , with τ = 0.5 , the autoimmune activation number is R 0 ( M 0 ) = 2.5287 .
The corresponding normalized sensitivity indices are summarized in Table 3 . 10.1371/journal.pone.0349131.t003 Table 3 Normalized sensitivity indices of the autoimmune activation number R 0 ( M 0 ) evaluated at the baseline parameter values. Parameter p Υ p R 0 η max 1.0000 θ − 0.3333 τ − 0.6667 M max 0.0564 α − 0.0047 δ 0.0047 ω − 0.0517 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 R 0 ( M 0 ) . In contrast, M max , ω , α , 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 R 0 ( M 0 ) 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 p = ( η max , θ , τ , M max , α , δ , ω ) , where each component is independently sampled within a prescribed interval around its baseline value p i * . Because experimentally established uncertainty ranges are unavailable for several parameters, a variation of ± 10 % is adopted as an exploratory uncertainty interval: p i ~ U ( 0.9 p i * , 1.1 p i * ) , i = 1 , . . . , 7 , where U denotes the uniform distribution.
A total of N s = 2000 parameter combinations are generated using LHS. For each sample, the autoimmune activation number R 0 ( M 0 ) is calculated using its analytical expression. The PRCC between a parameter p i and R 0 ( M 0 ) is computed as PRCC ( p i , R 0 ( M 0 ) ) = Corr ( r p i , r R 0 ) , where r p i and r R 0 denote the residuals obtained by regressing the rank-transformed values of p i and R 0 ( M 0 ) , 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 R 0 ( M 0 ) , 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. 10.1371/journal.pone.0349131.g006 Fig 6 Global sensitivity analysis of the autoimmune activation number R 0 ( M 0 ) using Latin Hypercube Sampling and Partial Rank Correlation Coefficients. The seven model parameters are independently varied within ± 10 % of their baseline values.
The maximum immune activation rate η max exhibits a strong positive correlation with R 0 ( M 0 ) , 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 M max 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, ‖ u ( k + 1 ) − u ( k ) ‖ < ε 1 , 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. 10.1371/journal.pone.0349131.g007 Fig 7 Optimal-control results for the multiple sclerosis model.
(A) Autoreactive immune-cell dynamics with and without treatment. (B) Cumulative objective value for the untreated and optimally controlled cases. (C) Optimal treatment intensity u ( t ).
(D) Comparison between the optimal treatment-adjustment rate a ( t ) and the numerical derivative du / dt , illustrating consistency with the treatment dynamics u ˙ ( t ) = a ( t ) . 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 A t = ∫ 0 t f I ( t ) d t , 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 F = ( 1 − A t A 0 ) × 100 , where A t is the cumulative autoimmune burden under optimal treatment and A 0 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 A 0 = 64.22 , whereas under optimal treatment it is reduced to A t = 1.73 . 10.1371/journal.pone.0349131.t004 Table 4 Effect of optimal treatment on cumulative autoimmune burden, treatment efficiency, and objective value.
Scenario Cumulative burden A Efficiency (%) Objective J No treatment 64.22 0 128.48 Optimal treatment 1.73 97.31 3.5676 This gives F = ( 1 − 1.73 64.22 ) × 100 ≈ 97.31 % . 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 J 0 = 128.48 without treatment to J t = 3.5676 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 u ˙ ( t ) = a ( t ) . 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 R 0 determines the local stability of the immune-free equilibrium, the analysis also reveals distinct dynamical regimes on either side of the critical threshold R 0 = 1 . In particular, when R 0 > 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 R 0 < 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 R 0 = 1 , the reduced system possesses a continuum of boundary equilibria, with E b = ( λ / ρ , 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 R 0 > 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 ].