A Type of Hopf Bifurcation in Models with Holling-II Type Multiple Sclerosis

Manhui He

College of Mathematics and Statistics, Northwest Normal University, Lanzhou, China.

DOI: 10.4236/jamp.2025.131015                

Abstract

Multiple autoimmune diseases often exhibit a cyclic pattern of relapse and remission, with significant periods of loss of self-tolerance being interrupted by recurrent autoimmune events. In this article, we explore a specific type of terminally differentiated regulatory T cell $HLA-D{R}^{+}{T}_{Reg}$ cells, and their application in existing autoimmune disease models. We also conduct an in-depth study on a multiple sclerosis model. This model incorporates a Holling-II type functional response mechanism. The focus of the study is to analyze whether the equilibrium points of the system have local asymptotic stability and determine the conditions for the existence of Hopf bifurcation. Furthermore, the direction of Hopf bifurcations and the stability of its periodic solutions can be analyzed through normal form theory and center manifold theorem.

Keywords

Multiple Sclerosis Model, Holling-II Type Functional Response, Equilibrium Point, Stability, Hopf Bifurcation

1. Introduction

Currently, in the medical field, autoimmune diseases have gradually become a focus of research and societal attention due to their significant impact on the health of populations worldwide. The immune system of the body not only serves as a faithful guardian to maintain health, but also is able to identify, track, and efficiently eliminate foreign pathogens in order to ensure balance and coordination within the internal environment. When facing autoimmune diseases, the immune system seems to have lost its proper recognition ability, mistakenly treating normal tissues and organs inside the body as foreign threats, thus initiating unnecessary immune responses that lead to a series of complex and diverse pathological changes and clinical symptoms.

Multiple sclerosis is a common autoimmune disease of the central nervous system, and its complex pathogenesis mainly stems from abnormal autoimmunity reactions. Once activated in the patient’s body, self-reactive T cells mistakenly identify myelin antigens of the central nervous system as foreign invaders. Subsequently, these T cells are able to penetrate the blood-brain barrier, move into the central nervous system, and launch an immune response to attack. This process results in the destruction of myelin sheath and damage to nerve fibers. B cells play a key role in this process, able to generate autologous antibodies specifically targeting myelin sheath components. Once these antibodies bind to myelin antigen, they can activate the complement system’s activity and exacerbate the damage of myelin sheath [1]. Multiple sclerosis involves a certain degree of genetic predisposition, with specific gene mutations closely related to its onset. For example, there is an association between specific alleles of the Human Leukocyte Antigen (HLA) gene, particularly HLA-DR15, and an increased likelihood of developing multiple sclerosis in individuals. The cell factor receptor genes such as IL-7RA have been confirmed to be associated with the occurrence of multiple sclerosis [2].

Recently, a novel effector cell subset $n{T}_{Reg}$ subgroup has been successfully identified in the research field. This unique subset of cells has the ability to further develop and differentiate into the final stage suppressor factors. Compared to traditional $n{T}_{Reg}$ cells, it performs better in terms of inhibitory function but has a relatively shorter lifespan. This discovery not only enriches our understanding of the immune regulation mechanism, but also provides a new direction for developing more precise and efficient immunotherapy strategies in the future. Phenotypic analysis reveals that the cell surface receptor HLA-DR encoded by specific genes exhibits diversity within T regulatory cells, successfully identifying and subclassifying different subgroups of terminal differentiation states [3]-[5]. Specifically, $HLA-D{R}^{+}T$ regulatory cells exhibit a faster ability to suppress the proliferation of conventional T cells compared to $HLA-D{R}^{-}T$ regulatory cells. According to our understanding, activation and expansion of $HLA-D{R}^{-}$ mediated effector nTreg cells drive the formation and development of this specific subset of $HLA-D{R}^{+}{T}_{Reg}$ cells.

Alexander and Wahl’s recent model successfully captured the internal feedback loop of autoimmunity. In this model, professional antigen-presenting cells ($pAPCs$ ) display self-antigens, triggering auto reactive effector T cells against host tissues and leading to attacks on the host tissue. Host tissue damage leads to an increased release of self-antigens, which in turn activates professional antigen-presenting cells ($pAPCs$ ). Tregs (regulatory T cells) play a crucial role in this process by suppressing autoreactive immune responses through multiple hypothesized mechanisms [6] [7].

For recent experimental results in the field of $HLA-D{R}^{+}{T}_{Reg}$ cells. We have decided to expand the model of W.J. ZHANG and LINDIM [8]. In the opening section, we presented a basic model and expanded upon it. In the second section, we obtained and analyzed the existence and stability of equilibrium points for the new model. In the third section. By selecting appropriate bifurcation parameters, we successfully identified the critical points of the Hopf bifurcation, which means that our new model can accurately.

2. Preliminaries

In 2014, W.J. ZHANG et al. [9] introduced regulatory T cells as an explicit variable into the model established by H.K. Alexander et al. obtaining the model:

$\{\begin{array}{l}\dot{A}=f\tilde{v}G-\left({\sigma }_{1}R+b_1\right)A-{\mu }_{A}A,\\ \dot{R}=\left({\pi }_{1}E+\beta \right)A-{\mu }_{R}R,\\ \dot{E}={\lambda }_{E}A-{\mu }_{E}E,\\ \dot{G}=\gamma E-\tilde{v}G-{\mu }_{G}G,\end{array}$ (2.1)

In this context, A, R, E and G respectively represent professional antigen-presenting cells (APCs), active regulatory T cells, activated effector T cells and specific self-antigens. $f\tilde{v}G$ represents the rate of maturation of professional APCs; ${\sigma }_{1}RA$ represents the rate at which active regulatory T cells inhibit professional APCs; activated auto reactive effector T are generated by resting T cells interacting with matured professional APCs at a rate of ${\lambda }_{E}A$ ; active regulatory T cell activation occurs at a rate $\left({\pi }_{1}E+\beta \right)A$ ; auto reactive effector T attack host tissues leading to release of self-antigens at a speed $\gamma E$ triggering maturation of professional $APCs$ starting a new cycle in autoimmunity. The death/clearance rates for populations A, R, E, and G are represented as ${\mu }_A$ , ${\mu }_R$ , ${\mu }_E$ and ${\mu }_G$ . Model (2.1) has an equilibrium point $E_0=\left(0,0,0,0\right)$ , and when $f\tilde{v}\gamma {\lambda }_E-\left(b_1+{\mu }_A\right){\mu }_E\left(\tilde{v}+{\mu }_G\right)>0$ , there is also has a stable equilibrium point $E_1=\left(A_1,R_1,E_1,G_1\right)$ , in which

$A_1=-\frac{\beta {\mu }_E}{2{\pi }_1{\lambda }_E}+\sqrt{{\left(\frac{\beta {\mu }_E}{2{\pi }_1{\lambda }_E}\right)}^2+\frac{f\tilde{v}\gamma {\lambda }_E-{\mu }_E\left(b_1+{\mu }_A\right)\left(\tilde{v}+{\mu }_G\right)}{{\sigma }_1{\pi }_1{\lambda }_E}}$

$R_1=\frac{f\tilde{v}\gamma {\lambda }_E-{\mu }_E\left(b_1+{\mu }_A\right)\left(\tilde{v}+{\mu }_G\right)}{{\sigma }_1{\pi }_1{\lambda }_E\left(\tilde{v}+{\mu }_G\right)},E_1=\frac{{\lambda }_E{A}_1}{{\mu }_E},G_1=\frac{\gamma {\lambda }_E{A}_1}{{\mu }_E\left(\tilde{v}+{\mu }_G\right)}.$

They provided the existence conditions and local stability of equilibrium points, but did not discuss the direction and stability of Hopf bifurcations. They only provide the conditions for the existence and local stability of equilibrium points, but do not discuss the direction and stability of Hopf bifurcations. Therefore, we introduce Holling-type functional response functions, which are important factors influencing system dynamics behavior, with Holling-II type functional response function being the most widely used. By studying systems with Holling-II type functional responses, analyzing the stability of all equilibrium points in the system and discovering that there is a Hopf bifurcation near a positive equilibrium point. Based on this foundation, we have obtained a new model.

$\{\begin{array}{l}\frac{\text{d}A}{\text{d}t}=f\tilde{v}G-\left(\frac{{\sigma }_{1}aR}{1+bR}+b_1\right)A-{\mu }_{A}A,\\ \frac{\text{d}R}{\text{d}t}=\left({\pi }_{1}E+\beta \right)A-{\mu }_{R}R,\\ \frac{\text{d}E}{\text{d}t}={\lambda }_{E}A-{\mu }_{E}E,\\ \frac{\text{d}G}{\text{d}t}=\gamma E-\tilde{v}G-{\mu }_{G}G,\end{array}$ (2.2)

Among them, $\frac{aR}{1+bR}$ is the Holling-II type functional response function. All

parameters in system (2.2) are positive numbers. For ease of calculation, we will perform a non-dimensionalization process on system (2.2), We order

$u=cA,v=dR,w=eE,x=yG,t=\frac{\tau }{q}.$

in which

$c=\frac{b\beta }{b_1+{\mu }_A},d=b,e=\frac{{\pi }_1}{\beta },y=\frac{\left(b_1+{\mu }_A\right){\pi }_1}{\beta \gamma },q=b_1+{\mu }_A.$

And still using $t$ instead of $\tau$ , the system (2.2) is reduced to

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=a_{1}x-\frac{a_{2}uv}{1+v}-u,\\ \frac{\text{d}v}{\text{d}t}=uw+u-a_{3}v,\\ \frac{\text{d}w}{\text{d}t}=b_{1}u-b_{2}w,\\ \frac{\text{d}x}{\text{d}t}=w-b_{3}x\end{array}$ (2.3)

in which

$a_1=\frac{b{\beta }^{2}f\tilde{v}\gamma }{{\pi }_1{q}^3},a_2=\frac{a{\sigma }_1}{bq},a_3=\frac{{\mu }_R}{q},b_1=\frac{{\pi }_1{\lambda }_E}{b{\beta }^2},b_2=\frac{{\mu }_E}{q},b_3=\frac{\tilde{v}+{\mu }_G}{q}.$

The parameters in system (2.3) are still positive. In terms of model parameter settings, each parameter has clear biological significance and rationale. For example, the parameter ${\sigma }_1$ is related to the ability of regulatory T cells to suppress professional APCs’ activity. Its magnitude reflects the inhibitory efficiency of regulatory T cells in the immune system, which may vary in different individuals or disease stages. Reasonable values are determined based on experimental research data on regulatory T cell function. Parameter ${\lambda }_E$ represents the rate at which resting T cells interact with mature professional APCs to generate activated effector T cells, and it is regulated by various cell factors and molecular signaling pathways. The value is determined based on experimental observations of relevant cell interactions. ${\mu }_A,{\mu }_R,{\mu }_E$ and ${\mu }_G$ represent the death/clearance rates of different cell populations, these rates are closely related to the lifespan, metabolic activity, and clearance mechanisms of the immune system in cells. Their reasonable ranges are determined through studying the survival cycles and metabolic processes of different cell types in vivo.

3. The Existence and Stability of Equilibrium Points

3.1. Existence of Equilibrium Point

For model (2.3),

(i) There always exists a trivial equilibrium point $E_0=\left(0,0,0,0\right)$ ;

(ii) only when

(H₁) $a_1{b}_1>b_2{b}_3$ and $a_1{b}_1-b_2{b}_3-a_2{b}_2{b}_3<0$ .

The system (2.3) has a unique positive equilibrium point. $E^{\text{*}}=\left(u^{\text{*}},v^{\text{*}},w^{\text{*}}{x}^{\text{*}}\right)$ , in which

$u^{\text{*}}=\frac{-A_2-\sqrt{A_2^2-4A_1{A}_3}}{2A_1},v^{\text{*}}=\frac{b_1{u}^{{\text{*}}^2}+b_2{u}^{\text{*}}}{a_3{b}_2},w^{\text{*}}=\frac{b_1{u}^{\text{*}}}{b_2},x^{\text{*}}=\frac{b_1{u}^{\text{*}}}{b_2{b}_3}.$

$A_1=a_1{b}_1^2-a_2{b}_1{b}_2{b}_3-b_1{b}_2{b}_3,A_2=a_1{b}_1{b}_2-a_2{b}_2^2{b}_3-b_3{b}_2^2,A_3=a_1{a}_3{b}_1{b}_2-a_3{b}_2^2{b}_3$

3.2. Local Stability of Equilibrium Point

The Jacobian matrix of system (2.3) at point $\left(u,v,w,x\right)$ is as follows

$J=\left(\begin{array}{cccc}-\frac{a_{2}v}{1+v}-1& -\frac{a_{2}u}{{\left(1+v\right)}^2}& 0& a_1\\ w+1& -a_3& u& 0\\ b_1& 0& -b_2& 0\\ 0& 0& 1& -b_3\end{array}\right),$ (3.1)

The stability of these equilibrium points is determined by calculating the eigenvalues of the Jacobian matrix at each equilibrium point using system (2.3).

Theorem 3.1. If condition $b_2{b}_3-a_1{b}_1>0$ , is satisfied, then the trivial equilibrium point $E_0=\left(0,0,0,0\right)$ is locally asymptotically stable; otherwise, $E_0$ is unstable.

Proof. The Jacobian matrix of system (2.3) at the equilibrium point $E_0$ is

$J_{E_0}=\left(\begin{array}{cccc}-1& 0& 0& a_1\\ 1& -a_3& 0& 0\\ b_1& 0& -b_2& 0\\ 0& 0& 1& -b_3\end{array}\right),$ (3.2)

The characteristic equation of matrix (3.2) is

${\lambda }^4+\left(h_1+1\right){\lambda }^3+\left(h_1+h_2\right){\lambda }^2+\left(h_2+h_3\right)\lambda +h_4.$ (3.3)

in which

$\begin{array}{l}{h}_1=a_3+b_2+b_3,h_2=a_3{b}_2+a_3{b}_3+b_2{b}_3,\\ h_3=a_3{b}_2{b}_3-a_1{b}_1,h_4=a_3{b}_2{b}_3-a_1{a}_3{b}_1\end{array}$

In order to ensure that all roots of the characteristic Equation (3.3) have negative real parts, we first calculate the following three sets of determinants [10].

${\text{Δ}}_1=h_1+1,{\text{Δ}}_2=\left|\begin{array}{cc}{h}_1+1& 1\\ h_2+h_3& h_1+h_2\end{array}\right|,{\text{Δ}}_3=\left|\begin{array}{ccc}{h}_1+1& 1& 0\\ h_2+h_3& h_1+h_2& h_1+1\\ 0& h_4& h_2+h_3\end{array}\right|$

Because ${\text{Δ}}_1,{\text{Δ}}_2$ is always greater than 0, to ensure that ${\text{Δ}}_3$ remains greater than 0, it is sufficient to satisfy condition $b_2{b}_3-a_1{b}_1>0$ . Therefore, when the conditions are met and judged by Routh-Hurwitz criteria, we can know that the roots of characteristic Equation (3.3) all have negative real parts, which means that the trivial equilibrium point $E_0$ of system (2.3) is locally asymptotically stable. Conversely, if not all roots of characteristic Equation (3.3) have negative real parts or if their real part is not zero, then the trivial equilibrium point $E_0$ will be unstable.

Theorem 3.2. Assume (H₁) holds. If satisfied

(H₂) $s_1,s_3>0$ , And ${\varphi }_3=\left(s_1{s}_2-s_3\right)s_3-s_1^2{s}_4>0$ .

The system (2.3) has a locally asymptotically stable equilibrium point $E^{\text{*}}$ , conversely, it is unstable when condition (H₂) is not satisfied.

Proof The Jacobian matrix of system (2.3) at the equilibrium point $E^{\text{*}}$ is

$J_{E^{\text{*}}}=\left(\begin{array}{cccc}-\frac{a_2{v}^{\text{*}}}{1+v^{\text{*}}}-1& -\frac{a_2{u}^{\text{*}}}{{\left(1+v^{\text{*}}\right)}^2}& 0& a_1\\ w^{\text{*}}+1& -a_3& u^{\text{*}}& 0\\ b_1& 0& -b_2& 0\\ 0& 0& 1& -b_3\end{array}\right),$ (3.4)

We order

$h_5=a_2,h_6=b_2+b_3,h_7=b_1$

$h_8=a_3{b}_2{b}_3,h_9=a_2{b}_3,h_{10}=a_2{b}_2{b}_3$

The characteristic equation of matrix (3.4) is

${\lambda }^4+s_1{\lambda }^3+s_2{\lambda }^2+s_3\lambda +s_4.$ (3.5)

$s_1=\frac{\left(h_1+h_5+1\right)v^{{\text{*}}^2}+\left(h_5+2h_1+2\right)v^{\text{*}}+\left(h_1+1\right)}{{\left(1+v^{\text{*}}\right)}^2},$

$s_2=\frac{\left(h_1{h}_5+h_1+h_2\right)v^{{\text{*}}^2}+\left(h_1{h}_5+2h_1+2h_2\right)v^{\text{*}}+h_1+h_2+h_1\left(u^{\text{*}}{w}^{\text{*}}+u^{\text{*}}\right)}{{\left(1+v^{\text{*}}\right)}^2},$

$s_3=\frac{\left(h_2{h}_5+h_2+h_3\right)v^{{\text{*}}^2}+\left(h_2{h}_5+2h_2+2h_3\right)v^{\text{*}}+h_2+h_3+h_5{h}_6\left(u^{\text{*}}{w}^{\text{*}}+u^{\text{*}}\right)}{{\left(1+v^{\text{*}}\right)}^2},$

$s_4=\frac{\left(h_5{h}_8+h_4\right)v^{{\text{*}}^2}+\left(h_5{h}_8+2h_4\right)v^{\text{*}}+h_4+h_{10}\left(u^{\text{*}}{w}^{\text{*}}+u^{\text{*}}\right)+h_7{h}_9{u}^{{\text{*}}^2}}{{\left(1+v^{\text{*}}\right)}^2}$

From this, it can be obtained

${\varphi }_1=s_1,{\varphi }_2=s_1{s}_2-s_3,{\varphi }_3=\left(s_1{s}_2-s_3\right)s_3-s_1^2{s}_4.$ (3.6)

Therefore, in the case where condition (H₂) is confirmed, according to the Routh-Hurwitz criterion, we can determine that all roots of characteristic Equation (3.5) have negative real parts, which means that the equilibrium point $E^{\text{*}}$ of system (2.3) exhibits local asymptotic stability. On the contrary, if not all roots of characteristic Equation (3.5) have negative real parts and their real parts are non-zero, then equilibrium point $E^{\text{*}}$ is unstable [11] [12].

Choose a few key parameters, such as $a_1,a_2$ and $a_3$ , analyze their effects on the stability of equilibrium points and system dynamics when they change within a certain range. Taking $a_1$ as an example, as $a_1$ increases gradually, it is observed that the stability of the trivial equilibrium point $E_0$ may change. The originally stable trivial equilibrium point may become unstable, and the stability of the positive equilibrium point $E^{\text{*}}$ will also be affected by this change with its stable region potentially shrinking or expanding. This is because the maturity rate of $a_1$ is related to professional APCs, and changes in $a_1$ will affect the dynamic balance of the entire immune system. Similarly, the changes in $a_2$ (related to the inhibitory effect of regulatory T cells on APC) will also lead to alterations in the stability of equilibrium points, affecting whether the system will exhibit a Hopf bifurcation and the parameter threshold for bifurcation [13]. Through this sensitivity analysis, we can gain a deeper understanding of the importance of various parameters in the model and how the system responds to changes in different biological processes.

4. Existence of Hopf Bifurcations

In this section, we select the parameter $h_5$ to study the existence of a Hopf branch of system (2.3) at the normal equilibrium point $E^{\text{*}}$ .

Theorem 4.1. Assuming that condition (H₁) holds, if there exists a $h_5^{*}>0$ , such that ${\varphi }_3\left(h_5^{*}\right)=\left(s_1\left(h_5^{*}\right)s_2\left(h_5^{*}\right)-s_3\left(h_5^{*}\right)\right)s_3\left(h_5^{*}\right)-s_1^2\left(h_5^{*}\right)s_4\left(h_5^{*}\right)=0$ and ${\frac{\text{d}\left(\varphi \left(h_5\right)\right)}{\text{d}{h}_5}|}_{h_5^{*}=h_5}\ne 0$ , then system (2.3) will exhibit a Hopf bifurcation at $h_5^{*}=h_5$ .

Proof. If there exists a $h_5^{\text{*}}>0$ such that ${\varphi }_3\left(h_5^{\text{*}}\right)=\left(s_1\left(h_5^{\text{*}}\right)s_2\left(h_5^{\text{*}}\right)-s_3\left(h_5^{\text{*}}\right)\right)s_3\left(h_5^{\text{*}}\right)-s_1^2\left(h_5^{\text{*}}\right)s_4\left(h_5^{\text{*}}\right)=0$ . Therefore, when $h_5^{\text{*}}=h_5$ , by characteristic Equation (3.5) can be written as

$\left({\lambda }^2+\frac{s_3}{s_1}\right)\left({\lambda }^2+s_1\lambda +s_2-\frac{s_3}{s_1}\right)=0.$ (4.1)

Suppose that $\lambda =\alpha \left(h_5\right)+\text{i}\beta \left(h_5\right)$ is a complex root of the characteristic Equation (2.5), substituting into Equation (3.5) we get

${\left(\alpha +\text{i}\beta \right)}^4+s_1{\left(\alpha +\text{i}\beta \right)}^3+s_2{\left(\alpha +\text{i}\beta \right)}^2+s_3\left(\alpha +\text{i}\beta \right)+s_4=0.$

When $h_5^{\text{*}}=h_5$ , $\alpha \left(h_5^{\text{*}}\right)=0$ . Taking the derivative of the above equation with respect to $h_5$ and simplifying, we get

$\begin{array}{l}-3{\alpha }^{\prime }\left(h_5^{\text{*}}\right)s_1\left(h_5^{\text{*}}\right){\beta }^2\left(h_5^{\text{*}}\right)-{s^{\prime }}_2\left(h_5^{\text{*}}\right){\beta }^2\left(h_5^{\text{*}}\right)-2s_2\left(h_5^{\text{*}}\right)\beta \left(h_5^{\text{*}}\right){\beta }^{\prime }\left(h_5^{\text{*}}\right)\\ +4{\beta }^3\left(h_5^{\text{*}}\right){\beta }^{\prime }\left(h_5^{\text{*}}\right)+s_3\left(h_5^{\text{*}}\right){\alpha }^{\prime }\left(h_5^{\text{*}}\right)+{s^{\prime }}_4\left(h_5^{\text{*}}\right)-4{\alpha }^{\prime }\left(h_5^{\text{*}}\right){\beta }^3\left(h_5^{\text{*}}\right)\text{i}\\ -{s^{\prime }}_1\left(h_5^{\text{*}}\right){\beta }^3\left(h_5^{\text{*}}\right)\text{i}-3s_1\left(h_5^{\text{*}}\right){\beta }^2\left(h_5^{\text{*}}\right){\beta }^{\prime }\left(h_5^{\text{*}}\right)\text{i}+2s_2\left(h_5^{\text{*}}\right){\alpha }^{\prime }\left(h_5^{\text{*}}\right)\beta \left(h_5^{\text{*}}\right)\text{i}\\ +{s^{\prime }}_3\left(h_5^{\text{*}}\right)\beta \left(h_5^{\text{*}}\right)\text{i}+s_3\left(h_5^{\text{*}}\right){\beta }^{\prime }\left(h_5^{\text{*}}\right)\text{i}=0\end{array}$

Separating the real part and the imaginary part, we can obtain

$\begin{array}{l}-3{\alpha }^{\prime }\left(h_5^{\text{*}}\right)s_1\left(h_5^{\text{*}}\right){\beta }^2\left(h_5^{\text{*}}\right)-{s^{\prime }}_2\left(h_5^{\text{*}}\right){\beta }^2\left(h_5^{\text{*}}\right)-2s_2\left(h_5^{\text{*}}\right)\beta \left(h_5^{\text{*}}\right){\beta }^{\prime }\left(h_5^{\text{*}}\right)\\ +4{\beta }^3\left(h_5^{\text{*}}\right){\beta }^{\prime }\left(h_5^{\text{*}}\right)+s_3\left(h_5^{\text{*}}\right){\alpha }^{\prime }\left(h_5^{\text{*}}\right)+{s^{\prime }}_4\left(h_5^{\text{*}}\right)=0,\end{array}$

$\begin{array}{l}(-4{\alpha }^{\prime }\left(h_5^{\text{*}}\right){\beta }^3\left(h_5^{\text{*}}\right)-{s^{\prime }}_1\left(h_5^{\text{*}}\right){\beta }^3\left(h_5^{\text{*}}\right)-3s_1\left(h_5^{\text{*}}\right){\beta }^2\left(h_5^{\text{*}}\right){\beta }^{\prime }\left(h_5^{\text{*}}\right)\\ +2s_2\left(h_5^{\text{*}}\right){\alpha }^{\prime }\left(h_5^{\text{*}}\right)\beta \left(h_5^{\text{*}}\right)+{s^{\prime }}_3\left(h_5^{\text{*}}\right)\beta \left(h_5^{\text{*}}\right)+s_3\left(h_5^{\text{*}}\right){\beta }^{\prime }\left(h_5^{\text{*}}\right))\text{i}=0.\end{array}$

Eliminating ${\beta }^{\prime }\left(h_5^{\text{*}}\right)$ from the above equation, we can obtain

${\alpha }^{\prime }\left(h_5^{\text{*}}\right)=\frac{N}{M},$

$\begin{array}{c}M=s_1^6\left(h_5^{\text{*}}\right)s_3\left(h_5^{\text{*}}\right)-6s_1^5\left(h_5^{\text{*}}\right)s_3^2\left(h_5^{\text{*}}\right)+4s_1^4\left(h_5^{\text{*}}\right)s_2^2\left(h_5^{\text{*}}\right)s_3\left(h_5^{\text{*}}\right)\\ +9s_1^4\left(h_5^{\text{*}}\right)s_3^3\left(h_5^{\text{*}}\right)-16s_1^4\left(h_5^{\text{*}}\right)s_2\left(h_5^{\text{*}}\right)s_3^3\left(h_5^{\text{*}}\right)+16s_3^5\left(h_5^{\text{*}}\right)\end{array}$

$\begin{array}{c}N=-s_1^6\left(h_5^{\text{*}}\right){s^{\prime }}_4\left(h_5^{\text{*}}\right)+3s_1^5\left(h_5^{\text{*}}\right)s_3\left(h_5^{\text{*}}\right){s^{\prime }}_4\left(h_5^{\text{*}}\right)-2s_1^2\left(h_5^{\text{*}}\right)s_2\left(h_5^{\text{*}}\right)s_3^2\left(h_5^{\text{*}}\right)\\ +s_1^4\left(h_5^{\text{*}}\right){s^{\prime }}_2\left(h_5^{\text{*}}\right)s_3^2\left(h_5^{\text{*}}\right)+2s_1^2\left(h_5^{\text{*}}\right){s^{\prime }}_1\left(h_5^{\text{*}}\right)s_2\left(h_5^{\text{*}}\right)s_3^3\left(h_5^{\text{*}}\right)\\ -3s_1^3\left(h_5^{\text{*}}\right){s^{\prime }}_2\left(h_5^{\text{*}}\right)s_3^3\left(h_5^{\text{*}}\right)+4s_1^2\left(h_5^{\text{*}}\right)s_3^3\left(h_5^{\text{*}}\right){s^{\prime }}_3\left(h_5^{\text{*}}\right)-4{s^{\prime }}_3\left(h_5^{\text{*}}\right)s_3^5\left(h_5^{\text{*}}\right)\end{array}$

And from (3.6), we get

$\begin{array}{c}\frac{\text{d}{\varphi }_3}{\text{d}{h}_5}={s^{\prime }}_1\left(h_5\right)s_2\left(h_5\right)s_3\left(h_5\right)+s_1\left(h_5\right){s^{\prime }}_2\left(h_5\right)s_3\left(h_5\right)+s_1\left(h_5\right)s_2\left(h_5\right){s^{\prime }}_3\left(h_5\right)\\ -2s_3\left(h_5\right){s^{\prime }}_3\left(h_5\right)-2s_1\left(h_5\right){s^{\prime }}_1\left(h_5\right)s_4\left(h_5\right)-s_1^2\left(h_5\right){s^{\prime }}_4\left(h_5\right).\end{array}$

Definition 4.2. Combining the expression of Theorem (4.1), we can obtain ${\frac{\text{d}{\varphi }_3}{\text{d}{h}_5}|}_{h_5=h_5^{*}}\ne 0$ , Therefore ${\frac{\text{d}\alpha \left(h_5\right)}{\text{d}{h}_5}|}_{h_5=h_5^{*}}\ne 0$ . Cross-sectional conditions hold, the system (2.3) will generate a $Hopf$ bifurcation at the positive equilibrium point $E^{*}$ .

The occurrence of Hopf bifurcation in this model has important biological significance. From a biological perspective, Hopf bifurcation means that the system transitions from a stable equilibrium state to periodic oscillations under specific parameter conditions. This may correspond to the transition of multiple sclerosis from a relatively stable period to a relapse period, where the dynamic balance of the immune system is disrupted, leading to changes in cell interactions and regulatory mechanisms, resulting in periodic fluctuations in disease-related indicators such as levels of self-reactive T cells and self-antigens. This oscillatory behavior may be related to the imbalance of feedback regulation within the immune system, such as disruption of the balance between regulatory T cells and effector T cells, leading to a periodic variation in disease progression [14].

5. Direction and Stability of Hopf Bifurcation

Hopf bifurcation theory focuses on the process in which stable equilibrium points transition from a stable state to a periodic oscillation state as system parameters change. Specifically to the normal number equilibrium point $E^{\text{*}}$ , first need to clarify the system’s differential equation group and determine the conditions for this point to be an equilibrium point. To prove the existence of a Hopf bifurcation at $E^{\text{*}}$ , it is typically necessary for the linearized equations of the system near equilibrium point $E^{\text{*}}$ to have eigenvalues with imaginary parts, indicating local instability. Nonlinear terms need to have an effect at a certain parameter value, allowing the system to transition from stable equilibrium states to periodic behavior. There exists a specific value of one or more system parameters, such that when these parameters cross this value, the behavior of the system transitions from stable equilibrium to periodic oscillations. The specific proof process may involve calculating the Jacobian matrix of the system at equilibrium point $E^{\text{*}}$ , analyzing its eigenvalues, and using tools such as center manifold theorem, Lyapunov exponents to determine if a Hopf bifurcation point exists. Through detailed mathematical derivation and analysis, it can be proven that under specific parameter conditions, the equilibrium point $E^{\text{*}}$ of normal numbers indeed exhibits a Hopf bifurcation phenomenon, revealing the complex dynamic behavior of the system as parameters change. In this section, we use the normal form theorem and projection method to discuss the direction and stability of the Hopf bifurcation generated near a non-zero equilibrium point $E^{\text{*}}$ of system (2.3) [15]. We first make a variable substitution for system (2.3).

$\tilde{u}=u-u^{\text{*}},\tilde{v}=v-v^{\text{*}},\tilde{w}=w-w^{\text{*}},\tilde{x}=x-x^{\text{*}}$

After the transformation, we still use $u,v,w,x$ instead of $\tilde{u},\tilde{v},\tilde{w},\tilde{x}$ , then system (2.3) becomes

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=a_1\left(x+x^{\text{*}}\right)-\frac{a_2\left(u+u^{\text{*}}\right)\left(v+v^{\text{*}}\right)}{1+\left(v+v^{\text{*}}\right)}-\left(u+u^{\text{*}}\right),\\ \frac{\text{d}v}{\text{d}t}=\left(u+u^{\text{*}}\right)\left(w+w^{\text{*}}\right)+\left(u+u^{\text{*}}\right)-a_3\left(v+v^{\text{*}}\right),\\ \frac{\text{d}w}{\text{d}t}=b_1\left(u+u^{\text{*}}\right)-b_2\left(w+w^{\text{*}}\right),\\ \frac{\text{d}x}{\text{d}t}=\left(w+w^{\text{*}}\right)-b_3\left(x+x^{\text{*}}\right).\end{array}$ (5.1)

Therefore, the positive equilibrium point $E^{\text{*}}=\left(u^{\text{*}},v^{\text{*}},w^{\text{*}},x^{\text{*}}\right)$ of system (2.3) becomes the zero point (0, 0, 0, 0) of system (4.1), then system (4.1) becomes

$\frac{\text{d}U}{\text{d}t}=J_{E^{*}}U+F\left(U\right)+\frac{1}{2}B\left(x,y\right)+\frac{1}{6}C\left(x,y,z\right)+\cdots ,$ (5.2)

Among them

$U={\left(u,v,w,x\right)}^{\text{T}},F\left(U\right)={\left(F_1,F_2,F_3,F_4\right)}^{\text{T}}.$

The components of $F\left(U\right)$ are as follows:

$F_1=-\frac{a_2}{{\left(1+v^{\text{*}}\right)}^2}uv+\frac{a_2{u}^{\text{*}}}{{\left(1+v^{\text{*}}\right)}^3}{v}^2+\frac{a_2}{{\left(1+v^{\text{*}}\right)}^3}u{v}^2-\frac{a_2{u}^{\text{*}}}{{\left(1+v^{\text{*}}\right)}^4}{v}^3+\cdots ,$

$F_2=uw+\cdots ,$

$F_3=0,$

$F_4=0.$

When $h_5=h_5^{\text{*}}$ , the eigenvalues of system (2.3) at equilibrium point $E^{\text{*}}$ are

${\lambda }_{1,2}=\frac{-s_1\pm \sqrt{s_1^2-4\left(s_2-\frac{s_3}{s_1}\right)}}{2},{\lambda }_{3,4}=\pm \text{i}{\omega }_0\left({\omega }_0=\sqrt{\frac{s_3}{s_1}}\right).$

Let ${\lambda }_{1,2}\ne 0$ , the eigenvalues of $J_E^{*}$ and $J_E^{*\text{T}}$ are purely imaginary roots corresponding eigenvectors as $q$ and $p$ meet

$J_E^{*}q=\text{i}{\omega }_{0}q,J_E^{*\text{T}}p=-\text{i}{\omega }_{0}p.$

$p=\left(\begin{array}{c}\frac{{\left(1+v^{\text{*}}\right)}^2\left[-a_3{b}_2+{\omega }_0^2+\left(a_3+b_2\right)\text{i}{\omega }_0\right]}{a_2{u}^{{\text{*}}^2}}\\ \frac{b_2-\text{i}{\omega }_0}{u^{\text{*}}}\\ 1\\ \frac{a_1{\left(1+v^{\text{*}}\right)}^2\left[-a_3{b}_2+{\omega }_0^2+\left(a_3+b_2\right)\text{i}{\omega }_0\right]}{a_2{u}^{{\text{*}}^2}\left(b_3-\text{i}{\omega }_0\right)}\end{array}\right),q=\left(\begin{array}{c}\frac{b_2{b}_3-{\omega }_0^2+\left(b_2+b_3\right)\text{i}{\omega }_0}{b_1}\\ \frac{b_2{b}_3+b_1{b}_3{u}^{\text{*}}-{\omega }_0^2+\left(b_2+b_3+b_1{u}^{\text{*}}\right)\text{i}{\omega }_0}{b_1\left(a_3+\text{i}{\omega }_0\right)}\\ b_3+\text{i}{\omega }_0\\ 1\end{array}\right).$

Let $q^{\text{*}}=q,\mu =\langle p,q\rangle ,p^{\text{*}}=p\frac{1}{\mu }$ . At this moment, satisfaction is achieved

$\mu ={\mu }_1+{\mu }_2+{\mu }_3+{\mu }_4,\langle p^{\text{*}},q^{\text{*}}\rangle =1,\langle {\overline{p}}^{\text{*}},q^{\text{*}}\rangle =\langle p^{\text{*}},{\overline{q}}^{\text{*}}\rangle =0.$

Among them

$\begin{array}{c}{\mu }_1=b_2^2{b}_3-a_3{b}_2{b}_3+a_3{\omega }_0^2+a_3{b}_2{\omega }_0^2+a_3{b}_3{\omega }_0^2+b_2^2{\omega }_0^2+2b_2{b}_3{\omega }_0^2-{\omega }_0^4\\ +\left(b_2^2+b_2{b}_3-a_3{b}_2-a_3{b}_3-a_3{b}_2{b}_3+a_3{\omega }_0^2-b_2^2{b}_3+2b_2{\omega }_0^2+b_3{\omega }_0^2\right)\text{i}{\omega }_0/a_2{b}_1{u}^{\ast 2},\end{array}$

$\begin{array}{c}{\mu }_2=b_2^2{b}_3+b_1{b}_2{b}_3{u}^{\text{*}}-2b_2{\omega }_0^2-b_3{\omega }_0^2-b_1{u}^{\text{*}}{\omega }_0^2\\ +\left(b_2^2+2b_2{b}_3+b_1{b}_2{u}^{\text{*}}+b_1{b}_2{b}_3{u}^{\text{*}}-{\omega }_0^2\right)\text{i}{\omega }_0/b_1{u}^{\text{*}}\left(a_3+\text{i}{\omega }_0\right),\end{array}$

${\mu }_3=b_3+\text{i}{\omega }_0,$

$\begin{array}{c}{\mu }_4={\left(1+v^{\text{*}}\right)}^2[-a_1{a}_3{b}_2{b}_3+\left(a_1{b}_3-a_1{a}_3-a_1{b}_2\right){\omega }_0^2\\ -\left(a_1{a}_3{b}_3+a_1{b}_2{b}_3-a_1{a}_3{b}_2+a_1{\omega }_0^2\right)\text{i}{\omega }_0]/a_2{u}^{\ast 2}\left(b_3^2+{\omega }_0^2\right).\end{array}$

${\overline{q}}^{*}=\left(\begin{array}{c}\frac{b_2{b}_3-{\omega }_0^2-\left(b_2+b_3\right)\text{i}{\omega }_0}{b_1}\\ \frac{a_3{b}_2{b}_3+a_3{b}_1{b}_3{u}^{\text{*}}+\left(b_1{u}^{\text{*}}+b_2+b_3-a_3\right){\omega }_0^2-\left(a_3{b}_2+a_3{b}_3+a_3{b}_1{u}^{\text{*}}-b_2{b}_3-b_1{b}_3+{\omega }_0^2\right)\text{i}{\omega }_0}{b_1\left(a_3^2+{\omega }_0^2\right)}\\ b_3-\text{i}{\omega }_0\\ 1\end{array}\right)$

${\overline{p}}^{\text{*}}$ and ${\overline{q}}^{\text{*}}$ are respectively the conjugate matrices of $p^{\text{*}}$ and $q^{\text{*}}$ , $\langle \rangle$ is the inner product symbol. According to the nonlinear term in system (5.2), consider a symmetric multilinear vector function with respect to $x={\left(x_1,x_2,x_3,x_4\right)}^{\text{T}}$ , $y={\left(y_1,y_2,y_3,y_4\right)}^{\text{T}}$ , $z={\left(z_1,z_2,z_3,z_4\right)}^{\text{T}}$ .

$B\left(x,y\right)=\left(\begin{array}{c}-\frac{a_2}{{\left(1+v^{\text{*}}\right)}^2}\left(x_1{y}_2+x_2{y}_1\right)+\frac{2a_2{u}^{\text{*}}}{{\left(1+v^{\text{*}}\right)}^3}{x}_2{y}_2\\ x_1{y}_3+x_3{y}_1\\ 0\\ 0\end{array}\right),$ (5.3)

$C\left(x,y,z\right)=\left(\begin{array}{c}\frac{2a_2{u}^{\text{*}}}{{\left(1+v^{\text{*}}\right)}^3}\left(x_1{y}_1{z}_2+x_1{y}_2{z}_1+x_2{y}_1{z}_1\right)\\ 0\\ 0\\ 0\end{array}\right).$ (5.4)

Definition 5.1. Decompose $U\in R^n$ into any component:

$U=z{q}^{\text{*}}+\overline{z}{\overline{q}}^{*}+y.$ (5.5)

The following relationship is given.

$\{\begin{array}{l}z=\langle p^{*},U\rangle ,\\ y=U-\langle p^{*},U\rangle q^{*}-\langle {\overline{p}}^{*},U\rangle {\overline{q}}^{*}\end{array}$ (5.6)

In the $\left(z,y\right)$ coordinates, Equation (5.2) is

$\{\begin{array}{l}\dot{z}=\text{i}{\omega }_{0}z+\langle p^{*},F\left(z{q}^{*}+\overline{z}{\overline{q}}^{*}+y\right)\rangle ,\\ \dot{y}=Ay+F\left(z{q}^{*}+\overline{z}{\overline{q}}^{*}+y\right)-\langle p^{*},F\left(z{q}^{*}+\overline{z}{\overline{q}}^{*}+y\right)\rangle q^{*}\\ -\langle {\overline{p}}^{*},F\left(z{q}^{*}+\overline{z}{\overline{q}}^{*}+y\right)\rangle {\overline{q}}^{*}\end{array}$ (5.7)

Using the properties of functions $B$ and $C$ , Equation (5.7) can be written as

$\{\begin{array}{l}\dot{z}=\text{i}{\omega }_{0}z+\frac{1}{2}{G}_{20}{z}^2+G_{11}z\overline{z}+\frac{1}{2}{G}_{02}{\overline{z}}^2+\frac{1}{2}{G}_{21}{z}^2\overline{z},\\ +\langle p^{*},B\left(q^{*},y\right)\rangle z+\langle p^{*},B\left({\overline{q}}^{*},y\right)\rangle \overline{z}+\cdots ,\\ \dot{y}=Ay+\frac{1}{2}{H}_{20}{z}^2+H_{11}z\overline{z}+\frac{1}{2}{H}_{02}{\overline{z}}^2+\cdots ,\end{array}$ (5.8)

Among them

$\begin{array}{l}{G}_{20}=\langle p^{*},B\left(q^{*},q^{*}\right)\rangle ,G_{11}=\langle p^{*},B\left(q^{*},{\overline{q}}^{*}\right)\rangle ,\\ G_{02}=\langle p^{*},B\left({\overline{q}}^{*},{\overline{q}}^{*}\right)\rangle ,G_{21}=\langle p^{*},C\left(q^{*},q^{*},{\overline{q}}^{*}\right)\rangle .\end{array}$ (5.9)

$\{\begin{array}{l}{H}_{20}=B\left(q^{*},q^{*}\right)-\langle p^{*},B\left(q^{*},q^{*}\right)\rangle q^{*}-\langle {\overline{q}}^{*},B\left(q^{*},q^{*}\right)\rangle {\overline{q}}^{*},\\ H_{11}=B\left(q^{*},{\overline{q}}^{*}\right)-\langle p^{*},B\left(q^{*},{\overline{q}}^{*}\right)\rangle q^{*}-\langle {\overline{q}}^{*},B\left(q^{*},{\overline{q}}^{*}\right)\rangle {\overline{q}}^{*},\\ H_{02}=B\left({\overline{q}}^{*},{\overline{q}}^{*}\right)-\langle p^{*},B\left({\overline{q}}^{*},{\overline{q}}^{*}\right)\rangle q^{*}-\langle {\overline{q}}^{*},B\left({\overline{q}}^{*},{\overline{q}}^{*}\right)\rangle {\overline{q}}^{*}\end{array}$ (5.10)

The calculation formula for the center manifold in Equation (5.8) is

$w_{20}={\left(2\text{i}{\omega }_{0}I-A\right)}^{-1}{H}_{20},w_{11}=-A^{-1}{H}_{11},w_{02}={\left(-2\text{i}{\omega }_{0}I-A\right)}^{-1}{H}_{02}.$ (5.11)

where $I$ is the n-order identity matrix, substituting the central manifold $y=w\left(z,\overline{z}\right)$ into the first equation of formula (5.8), We obtain the following simplified equation.

$\begin{array}{c}\dot{z}=\text{i}{\omega }_{0}z+\frac{1}{2}{G}_{20}{z}^2+G_{11}z\overline{z}+\frac{1}{2}{G}_{02}{\overline{z}}^2+\frac{1}{2}(G_{21}-2\langle p^{*},B\left(q^{*},A^{-1}{H}_{11}\right)\rangle \\ +\langle p^{*},B\left({\overline{q}}^{*},{\left(2\text{i}{\omega }_{0}I-A\right)}^{-1}{H}_{20}\right)\rangle z^2\overline{z}+\cdots \end{array}$ (5.12)

or

$\dot{z}=\text{i}{\omega }_{0}z+\frac{1}{2}{g}_{20}{z}^2+g_{11}z\overline{z}+\frac{1}{2}{g}_{02}{\overline{z}}^2+\frac{1}{2}{g}_{21}{z}^2\overline{z}+\cdots ,$ (5.13)

Among them $g_{20}=G_{20},g_{11}=G_{11},g_{02}=G_{02}$ . And

$g_{21}=G_{21}-2\langle p^{*},B\left(q^{*},A^{-1}{H}_{11}\right)\rangle +\langle p^{*},B\left({\overline{q}}^{*},{\left(2\text{i}{\omega }_{0}I-A\right)}^{-1}{H}_{20}\right)\rangle .$

Using the following relationship

$A^{-1}{q}^{*}=\frac{q^{*}}{\text{i}{\omega }_0},A^{-1}{\overline{q}}^{*}=-\frac{{\overline{q}}^{*}}{\text{i}{\omega }_0},{\left(2\text{i}{\omega }_{0}I-A\right)}^{-1}{q}^{*}=\frac{q^{*}}{i{\omega }_0},{\left(2\text{i}{\omega }_{0}I-A\right)}^{-1}{\overline{q}}^{*}=\frac{{\overline{q}}^{*}}{3\text{i}{\omega }_0},$

Simplify A to

$\begin{array}{c}{g}_{21}=\langle p^{*}C\left(q^{*},q^{*},{\overline{q}}^{*}\right)\rangle -2\langle p^{*},B\left(q^{*},A^{-1}B\left(q^{*},{\overline{q}}^{*}\right)\right)\rangle \\ +\langle p^{*},B\left({\overline{q}}^{*},{\left(2\text{i}{\omega }_{0}I-A\right)}^{-1}B\left(q^{*},q^{*}\right)\right)\rangle \\ +\frac{1}{\text{i}{\omega }_0}\langle p^{*},B\left(q^{*},q^{*}\right)\rangle \langle p^{*},B\left(q^{*},{\overline{q}}^{*}\right)\rangle \\ -\frac{2}{\text{i}{\omega }_0}{\left|\langle p^{*},B\left(q^{*},{\overline{q}}^{*}\right)\rangle \right|}^2-\frac{1}{3\text{i}{\omega }_0}{\left|\langle p^{*},B\left({\overline{q}}^{*},{\overline{q}}^{*}\right)\rangle \right|}^2.\end{array}$ (5.14)

The formula for calculating the first Lyapunov coefficient is as follows:

$\begin{array}{c}{l}_1\left(0\right)=\frac{1}{2{\omega }_0}\text{Re}\{\langle p^{\text{*}},C\left(q^{\text{*}},q^{\text{*}},{\overline{q}}^{\text{*}}\right)\rangle -2\langle p^{\text{*}},B\left(q^{\text{*}},A^{-1}B\left(q^{\text{*}},{\overline{q}}^{\text{*}}\right)\right)\rangle \\ +\langle p^{\text{*}}B\left({\overline{q}}^{\text{*}},{\left(2\text{i}{\omega }_0{I}_4-A\right)}^{-1}B\left(q^{\text{*}},q^{\text{*}}\right)\right)\rangle \}\end{array}$ (5.15)

Among them $A=J_E^{*}$ .

According to the expression of (5.5), we can calculate

$B\left(q^{\text{*}},q^{\text{*}}\right)=\left(\begin{array}{c}-\frac{2a_2}{{\left(1+v^{\text{*}}\right)}^2}{q}_1{q}_2+\frac{2a_2{u}^{\text{*}}}{{\left(1+v^{\text{*}}\right)}^3}{q}_2{}^2\\ 2q_1{q}_3\\ 0\\ 0\end{array}\right),$ (5.16)

$B\left(q^{\text{*}},{\overline{q}}^{*}\right)=\left(\begin{array}{c}-\frac{a_2}{{\left(1+v^{\text{*}}\right)}^2}\left(q_1{\overline{q}}_2+q_2{\overline{q}}_1\right)+\frac{2a_2{u}^{\text{*}}}{{\left(1+v^{\text{*}}\right)}^3}{q}_2{\overline{q}}_2\\ q_1{\overline{q}}_3+q_3{\overline{q}}_1\\ 0\\ 0\end{array}\right),$ (5.17)

$C\left(q^{\text{*}},q^{\text{*}},{\overline{q}}^{*}\right)=\left(\begin{array}{c}\frac{2a_2{u}^{\text{*}}}{{\left(1+v^{\text{*}}\right)}^3}\left(q_1{}^2{\overline{q}}_2+2q_1{q}_2{\overline{q}}_1\right)\\ 0\\ 0\\ 0\end{array}\right),$ (5.18)

$q_1=\frac{b_2{b}_3-{\omega }_0^2+\left(b_2+b_3\right)\text{i}{\omega }_0}{b_1},$

$q_2=\frac{b_2{b}_3+b_1{b}_3{u}^{\text{*}}-{\omega }_0^2+\left(b_2+b_3+b_1{u}^{\text{*}}\right)\text{i}{\omega }_0}{b_1\left(a_3+\text{i}{\omega }_0\right)},$

$q_3=b_3+\text{i}{\omega }_0,$

${\overline{q}}_1=\frac{b_2{b}_3-{\omega }_0^2-\left(b_2+b_3\right)\text{i}{\omega }_0}{b_1},$

${\overline{q}}_2=\frac{a_3{b}_2{b}_3+a_3{b}_1{b}_3{u}^{\text{*}}+\left(b_1{u}^{\text{*}}+b_2+b_3-a_3\right){\omega }_0^2-\left(a_3{b}_2+a_3{b}_3+a_3{b}_1{u}^{\text{*}}-b_2{b}_3-b_1{b}_3+{\omega }_0^2\right)i{\omega }_0}{b_1\left(a_3^2+{\omega }_0^2\right)},$

${\overline{q}}_3=b_3-\text{i}{\omega }_0.$

$A^{-1}=\frac{1}{\kappa }\left|\begin{array}{llll}{A}_{11}\hfill & A_{12}\hfill & A_{13}\hfill & A_{14}\hfill \\ A_{21}\hfill & A_{22}\hfill & A_{23}\hfill & A_{24}\hfill \\ A_{31}\hfill & A_{32}\hfill & A_{33}\hfill & A_{34}\hfill \\ A_{41}\hfill & A_{42}\hfill & A_{43}\hfill & A_{44}\hfill \end{array}\right|,{\left(2i{\omega }_0{I}_4-A\right)}^{-1}=\frac{1}{\psi }\left|\begin{array}{llll}{B}_{11}\hfill & B_{12}\hfill & B_{13}\hfill & B_{14}\hfill \\ B_{21}\hfill & B_{22}\hfill & B_{23}\hfill & B_{24}\hfill \\ B_{31}\hfill & B_{32}\hfill & B_{33}\hfill & B_{34}\hfill \\ B_{41}\hfill & B_{42}\hfill & B_{43}\hfill & B_{44}\hfill \end{array}\right|$

$\kappa =u^{\text{*}}{b}_1{b}_3{d}_2+a_3{b}_2{b}_3{d}_1+b_2{b}_3{d}_2{d}_3+a_1{a}_3{b}_1.$

$\begin{array}{c}\psi =-16{\omega }_0^4+\left(8\text{i}{a}_3+8\text{i}{b}_2+8\text{i}{b}_3-8\text{I}{d}_1\right){\omega }_0^3+(4a_3{b}_2+4a_3{b}_3-4a_3{d}_1+4b_2{b}_3\\ -4b_2{d}_1-4b_3{d}_1-4d_2{d}_3){\omega }_0^2+(2\text{i}{u}^{\text{*}}{b}_1{d}_2-2\text{i}{a}_3{b}_2{b}_3+2\text{i}{a}_3{b}_2{d}_1+2\text{i}{a}_3{b}_3{d}_1\\ +2\text{i}{b}_2{b}_3{d}_1+2\text{i}{b}_2{d}_2{d}_3+2\text{i}{b}_3{d}_2{d}_3+2\text{i}{a}_1{b}_1){\omega }_0+u^{\text{*}}{b}_1{b}_3{d}_2+a_3{b}_2{b}_3{d}_1\\ +d_3{d}_2{b}_2{b}_3+a_1{b}_1{a}_3.\end{array}$

$A_{11}=a_3{b}_2{b}_3,A_{12}=b_2{b}_3{d}_2,A_{13}=u^{\text{*}}{b}_3{d}_2+a_1{a}_3,A_{14}=a_1{a}_3{b}_2,$

$A_{21}=b_3\left(u^{\text{*}}{b}_1+b_2{d}_3\right),A_{22}=-\left(b_2{b}_3{d}_1+a_1{b}_1\right),A_{23}=-\left(u^{\text{*}}{b}_3{d}_1-a_1{d}_3\right),$

$A_{24}=a_1\left(u^{\text{*}}{b}_1+b_2{d}_3\right),A_{31}=a_3{b}_1{b}_3,A_{32}=b_1{b}_3{d}_2,$

$A_{33}=-\left(a_3{b}_3{d}_1+a_3{d}_2{d}_3\right),A_{34}=a_1{a}_3{b}_1,A_{41}=a_3{b}_1,A_{42}=b_1{d}_2,$

$A_{43}=-\left(a_3{d}_1+d_2{d}_3\right),A_{44}=-\left(u^{\text{*}}{b}_1{d}_2+a_3{b}_2{d}_1+b_2{d}_2{d}_3\right),$

$B_{11}=-\left(2\text{i}{\omega }_0+a_3\right)\left(2\text{i}{\omega }_0+b_2\right)\left(2\text{i}{\omega }_0+b_3\right),B_{12}=-d_2\left(2\text{i}{\omega }_0+b_2\right)\left(2\text{i}{\omega }_0+b_3\right),$

$B_{13}=-2\text{i}{u}^{\text{*}}{d}_2{\omega }_0+2\text{i}{a}_1{\omega }_0+u^{\text{*}}{b}_3{d}_2+a_1{a}_3,B_{14}=-a_1\left(2\text{i}{\omega }_0+a_3\right)\left(2\text{i}{\omega }_0+b_2\right),$

$B_{21}=4d_3{\omega }_0^2-2\text{i}{u}^{\text{*}}{b}_1{\omega }_0-2\text{i}{b}_2{d}_3{\omega }_0-2\text{i}{b}_3{d}_3{\omega }_0-u^{\text{*}}{b}_1{b}_3-b_2{b}_3{d}_3,$

$\begin{array}{c}{B}_{22}=8\text{i}{\omega }_0^3+4b_2{\omega }_0^2+4b_3{\omega }_0^2-4d_1{\omega }_0^2-2\text{i}{b}_2{b}_3{\omega }_0\\ +2\text{i}{b}_2{d}_1{\omega }_0+2\text{i}{b}_3{b}_1{\omega }_0+b_2{b}_3{d}_1+a_1{b}_1,\end{array}$

$B_{23}=4u^{\text{*}}{\omega }_0^2-2\text{i}{u}^{\text{*}}{b}_3{\omega }_0+2\text{i}{u}^{\text{*}}{d}_1{\omega }_0+u^{\text{*}}{b}_3{d}_1-d_3{a}_1,$

$B_{24}=-a_1\left(2\text{i}{u}^{\text{*}}{d}_3{\omega }_0+u^{\text{*}}{b}_1+b_2{d}_3\right),$

$B_{31}=-b_1\left(2\text{i}{\omega }_0+a_3\right)\left(2\text{i}{\omega }_0+b_3\right),B_{32}=-b_1{d}_2\left(2\text{i}{\omega }_0+b_3\right),$

$\begin{array}{c}{B}_{33}=8\text{i}{\omega }_0^3+4a_3{\omega }_0^2+4b_3{\omega }_0^2-4d_1{\omega }_0^2-2\text{i}{a}_3{b}_3{\omega }_0+2\text{i}{a}_3{d}_1{\omega }_0\\ +2\text{i}{b}_3{d}_1{\omega }_0+2\text{i}{d}_2{d}_3{\omega }_0+a_3{b}_3{d}_1+b_3{d}_2{d}_3,\end{array}$

$B_{34}=-a_1{b}_1\left(2\text{i}{\omega }_0+a_3\right),B_{41}=-b_1\left(2\text{i}{\omega }_0+a_3\right),$

$B_{42}=-b_1{d}_2,B_{43}=4{\omega }_0^2-2\text{i}{a}_3{\omega }_0+2\text{i}{d}_1{\omega }_0+a_3{d}_1-d_2{d}_3,$

$\begin{array}{c}{B}_{44}=8\text{i}{\omega }_0^3+4a_3{\omega }_0^2+4b_2{\omega }_0^2-4d_1{\omega }_0^2-2\text{i}{a}_3{b}_2{\omega }_0+2\text{i}{a}_3{d}_1{\omega }_0\\ +2\text{i}{b}_2{d}_1{\omega }_0+2\text{i}{d}_2{d}_3{\omega }_0+u^{\text{*}}{b}_1{d}_2+a_3{b}_2{d}_1+b_2{d}_2{d}_3,\end{array}$

$d_1=-\frac{a_2{v}^{\text{*}}}{1+v^{\text{*}}}-1,d_2=-\frac{a_2{u}^{\text{*}}}{{\left(1+v^{\text{*}}\right)}^2},d_3=w^{\text{*}}+1.$

Finally, from (5.15) we can obtain the expression for the first Lyapunov coefficient. Given that the expression of the first Lyapunov coefficient is complex and lengthy, to avoid redundancy, the specific form will not be detailed here. However, we will effectively verify the directionality and stability of the Hopf bifurcation through numerical simulation methods in the conclusion section [16].

Theorem 5.2. Assuming that (H₁)-(H₂) holds, the system (2.3) will generate a Hopf bifurcation at $h_5=h_5^{*}$ .

(i) If $l_1\left(0\right)<0$ , then the branch periodic solutions of Hopf bifurcation are orbitally asymptotically stable and branching direction is subcritical.

(ii) If $l_1\left(0\right)>0$ , then the branch periodic solutions of Hopf bifurcation are unstable and the branching direction is supercritical.

In the research of multiple sclerosis, there are many experiments and clinical data available to validate the effectiveness and relevance of the model. For example, it has been clinically observed that patients with multiple sclerosis show a significant increase in the number of self-reactive T cells in their blood and cerebrospinal fluid during disease relapse, along with elevated levels of myelin-specific antibodies, which is accompanied by worsening neurological dysfunction. In the model, when a Hopf bifurcation occurs leading to periodic oscillations, there is a corresponding fluctuation in the number of self-reactive T cells and levels of self-antigens, which matches with clinical observations of relapse period characteristics.

From experimental data, some studies have found that in animal models of multiple sclerosis, when the regulatory T cell function is impaired or decreased, the disease worsens and is more likely to relapse. In our model, changes in parameters affect the function of regulatory T cells (such as parameter $a_2$ being related to regulatory T cells). When these parameter changes result in stability alterations or Hopf bifurcations at equilibrium points, it reflects the impact of changes in regulatory T cell function on the entire immune system, consistent with experimental results. Through this comparison, it further proves that the model can simulate to a certain extent the pathophysiological process of multiple sclerosis, providing valuable theoretical basis for the research and treatment of the disease.

6. Numerical Simulation

Select parameters $a_1=2.5$ , $a_2=1.5$ , $a_3=0.006$ , $b_1=2$ , $b_2=1$ , $b_3=3.9$ , The equilibrium point at this moment is

$E^{\text{*}}=\left(0.001385633737,0.2315789497,0.002771267474,0.0007105814036\right).$

$a_1{b}_1-b_2{b}_3=1.1>0$ , $a_1{b}_1-b_2{b}_3-a_2{b}_2{b}_3=-4.75<0$ . $h_1=4.906$ , $h_2=3.9294$ , $h_3=-4.9766$ , $h_4=-0.0066$ , $s_1=6.188051281$ , $s_2=10.22363780$ , $s_3=0.06782538726$ , $s_4=0.005373784551$ . ${\text{Δ}}_1=5.906>0$ , ${\text{Δ}}_2=53.2290724>0$ , ${\text{Δ}}_3=-55.51127110<0$ . ${\varphi }_1=6.188051281>0$ , ${\varphi }_2=63.19656959>0$ , ${\varphi }_3=4.080558963>0$ . Therefore, the trivial equilibrium point of system (2.3) is unstable, while the positive equilibrium point is locally asymptotically stable (see Figure 1). Take parameter $a_2=0.289$ , $a_3=0.00001$ again, while keeping the rest of the parameter values unchanged, then the equilibrium point becomes

$E^{\text{*}}=\left(0.0004055750922,40.59040744,0.0008111501844,0.0002079872268\right).$

And $s_1=6.182061281$ , $s_2=10.18211425$ , $s_3=0.0001021538119$ , $s_4=2.646984371\times {10}^{-7}$ . ${\varphi }_1=6.182061281>0$ , ${\varphi }_2=62.94635211>0$ , ${\varphi }_3=\left(s_1{s}_2-s_3\right)s_3-s_1^2{s}_4\approx 0$ , ${\alpha }^{\prime }\left(h_5^{\text{*}}\right)=-55.6904867<0$ . Calculate obtained

$p^{*}=\left(\begin{array}{c}-1.301726+0.559858\text{i}\\ 9.385868\times {10}^{-6}+0.000021\text{i}\\ 3.770875\times {10}^{-9}+8.806580\times {10}^{-9}\text{i}\\ -0.834813+0.358013\text{i}\end{array}\right)\text{,}{q}^{*}=\left(\begin{array}{c}1.949991+0.009959\text{i}\\ 1.951573+0.009960\text{i}\\ 3.9+0.004065\text{i}\\ 1\end{array}\right),$

This can be obtained from (5.16)-(5.18).

$B\left(q^{*},q^{*}\right)=\left(\begin{array}{c}-0.001271-0.000012\text{i}\\ 15.209854+0.093535\text{i}\\ 0\\ 0\end{array}\right),B\left(q^{*},{\overline{q}}^{*}\right)=\left(\begin{array}{c}-0.000635-7.745295\times {10}^{-6}\text{i}\\ 15.210016\\ 0\\ 0\end{array}\right)$

$C\left(q^{*},q^{*},{\overline{q}}^{*}\right)=\left(\begin{array}{c}4.836052568\times {10}^{-8}+5.414904080\times {10}^{-10}\text{i}\\ 0\\ 0\\ 0\end{array}\right).$

Further obtain

$\begin{array}{l}\langle p^{*},C\left(q^{*},q^{*},{\overline{q}}^{*}\right)\rangle =-6.325532836\times {10}^{-8}+2.637015886\times {10}^{-8}\text{i},\\ 2\langle p^{*},B\left(q^{*},A^{-1}B\left(q^{*},{\overline{q}}^{*}\right)\right)\rangle =8.032371924-3.285379051\text{i},\\ \langle p^{*},B\left({\overline{q}}^{*},{\left(2\text{i}{\omega }_0{I}_4-A\right)}^{-1}B\left(q^{*},q^{*}\right)\right)\rangle =-113.9881448-174.6288817\text{i},\end{array}$

Further substituting the above results into (5.15) yields $l_1\left(0\right)=-13032.68784<0$ . by Theorem 4 it is known that a subcritical Hopf bifurcation occurs at the normal form equilibrium point $E^{\text{*}}$ and the periodic solutions on this branch are stable (see Figure 2). Select parameters again

$a_1=1.5,a_2=3,a_3=0.0000225,b_1=1,b_2=0.6,b_3=2.$

By using the same method as above, it can be calculated that $l_1\left(0\right)=8.740490825\times {10}^6>0$ . According to Theorem 4, there is a subcritical Hopf branch generated at the equilibrium point $E^{\text{*}}$ of the normal number (see Figure 3), and the periodic solution of the branch is unstable. Additionally, take appropriate parameter values

(a)

(b)

Figure 1. The graphical representation and trajectory of numerical solutions $u\left(t\right),v\left(t\right),w\left(t\right),x\left(t\right)$ when parameter $h_5=1.5$ .

(a)

(b)

Figure 2. The graphical representation and trajectory of numerical solutions $u\left(t\right),v\left(t\right),w\left(t\right),x\left(t\right)$ when parameter $h_5=0.289$ .

$a_1=0.15,a_2=0.3,a_3=0.002,b_1=0.1,b_2=0.06,b_3=\mathrm{0.2.}$

The system generates limit cycles (see Figure 4).

Where $u\left(t\right),v\left(t\right),w\left(t\right)$ and $x\left(t\right)$ represent the variables of antigen-presenting cells (Apcs), active regulatory T cells, active effector T cells, and specific self-antigens changing over time. The trajectory diagram of $u\left(t\right),v\left(t\right),w\left(t\right)$ and $x\left(t\right)$ within the range in three-dimensional graphics.

Figure 3. The graphical representation and trajectory of numerical solutions $u\left(t\right),v\left(t\right),w\left(t\right),x\left(t\right)$ when parameter $h_5=3$ .

(a)

(b)

Figure 4. The graph of numerical solutions $u\left(t\right),v\left(t\right),w\left(t\right),x\left(t\right)$ when parameter $h_5=0.3$ , and the system has limit cycles.

7. Conclusions

In this paper, we discuss a multiple sclerosis model with Holling-II type functional reactions. The system (2.3) always has a trivial equilibrium point $E_0\left(0,0,0,0\right)$ , and then we discuss the conditions for the existence of positive constant equilibrium points in the model, and study the stability of the system at trivial equilibrium points and positive equilibrium points. Proved that when $s_1,s_3>0$ is satisfied, and the system’s equilibrium point is locally asymptotically stable when also satisfying ${\varphi }_3=\left(s_1{s}_2-s_3\right)s_3-s_1^2{s}_4>0$ . Furthermore, by selecting the rate of inhibition cells $T$ as a parameter for regulating active cell signaling in antigen-presenting cells, when $h_5$ increases to a critical value $h_5^{\text{*}}$ , system (2.3) will generate a Hopf bifurcation at the steady state point $E^{\text{*}}$ .

In the context of multiple sclerosis, the negative feedback mechanism between effector $T$ cells and regulatory T cells is considered a potential checkpoint for regulation. However, our research revealed that increasing the activity of T cells or self-reactive T cells can cause damage, and upregulating $\alpha$ parameters significantly increases the risk of recurrence. Increase the parameter ${\lambda }_E$ , that is, enhancing the efficiency of T cell activation by antigen-presenting cells specialized in presenting antigens, also produces similar effects. This model emphasizes the importance of the positive feedback loop between enhanced professional antigen-presenting cells and self-reactive effector T cells, offering a promising path to be a key strategy in preventing relapses of multiple sclerosis. With a deeper understanding of the mechanism of regulatory T cells, it may be possible in the future to achieve more precise predictions of relapses in multiple sclerosis.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Alexander, H.K. and Wahl, L.M. (2010) Self-Tolerance and Autoimmunity in a Regulatory T Cell Model. Bulletin of Mathematical Biology, 73, 33-71. https://doi.org/10.1007/s11538-010-9519-2
[2]	Stürner, K.H., Siembab, I., Schön, G., Stellmann, J., Heidari, N., Fehse, B., et al. (2019) Is Multiple Sclerosis Progression Associated with the HLA-DR15 Haplotype? Multiple Sclerosis Journal—Experimental, Translational and Clinical, 5. https://doi.org/10.1177/2055217319894615
[3]	Caridade, M., Graca, L. and Ribeiro, R.M. (2013) Mechanisms Underlying CD4+ Treg Immune Regulation in the Adult: From Experiments to Models. Frontiers in Immunology, 4, Article No. 378. https://doi.org/10.3389/fimmu.2013.00378
[4]	Moser, T., Hoepner, L., Schwenker, K., Seiberl, M., Feige, J., Akgün, K., et al. (2021) Cladribine Alters Immune Cell Surface Molecules for Adhesion and Costimulation: Further Insights to the Mode of Action in Multiple Sclerosis. Cells, 10, Article No. 3116. https://doi.org/10.3390/cells10113116
[5]	Baecher‐Allan, C. and Hafler, D.A. (2006) Human Regulatory T Cells and Their Role in Autoimmune Disease. Immunological Reviews, 212, 203-216. https://doi.org/10.1111/j.0105-2896.2006.00417.x
[6]	Putnam, A.L., Brusko, T.M., Lee, M.R., Liu, W., Szot, G.L., Ghosh, T., et al. (2009) Expansion of Human Regulatory T-Cells from Patients with Type 1 Diabetes. Diabetes, 58, 652-662. https://doi.org/10.2337/db08-1168
[7]	Roncador, G., Brown, P.J., Maestre, L., Hue, S., Martínez-Torrecuadrada, J.L., Ling, K., et al. (2005) Analysis of FOXP3 Protein Expression in Human CD4+CD25+ Regulatory T Cells at the Single-Cell Level. European Journal of Immunology, 35, 1681-1691. https://doi.org/10.1002/eji.200526189
[8]	Zhang, W., Wahl, L.M. and Yu, P. (2014) Modeling and Analysis of Recurrent Autoimmune Disease. SIAM Journal on Applied Mathematics, 74, 1998-2025. https://doi.org/10.1137/140955823
[9]	Zhang, W., Wahl, L.M. and Yu, P. (2014) Viral Blips May Not Need a Trigger: How Transient Viremia Can Arise in Deterministic In-Host Models. SIAM Review, 56, 127-155. https://doi.org/10.1137/130937421
[10]	Yu, P. (2005) Closed-Form Conditions of Bifurcation Points for General Differential Equations. International Journal of Bifurcation and Chaos, 15, 1467-1483. https://doi.org/10.1142/s0218127405012582
[11]	Jeschke, J.M., Kopp, M. and Tollrian, R. (2002) Predator Functional Responses: Discriminating between Handling and Digesting Prey. Ecological Monographs, 72, 95-112. https://doi.org/10.1890/0012-9615(2002)072[0095:pfrdbh]2.0.co;2
[12]	Dawes, J.H.P. and Souza, M.O. (2013) A Derivation of Holling’s Type I, II and III Functional Responses in Predator-Prey Systems. Journal of Theoretical Biology, 327, 11-22. https://doi.org/10.1016/j.jtbi.2013.02.017
[13]	Hinrichsen, D. and Pritchard, A.J. (2005) Mathematical Systems Theory I: Modelling, State Space Analysis, Stability and Robustness. Springer.
[14]	Lucchinetti, C., Brück, W., Parisi, J., Scheithauer, B., Rodriguez, M. and Lassmann, H. (2000) Heterogeneity of Multiple Sclerosis Lesions: Implications for the Pathogenesis of Demyelination. Annals of Neurology, 47, 707-717. https://doi.org/10.1002/1531-8249(200006)47:6<707::aid-ana3>3.0.co;2-q
[15]	Kuznetsov, Y.A. (2023) Elements of Applied Bifurcation Theory. Springer.
[16]	Hassard, B.D., Kazarinoff, N.D. and Wan, Y.H. (1981) Theory and Applications of Hopf Bifurcation. Cambridge University Press.