Nonlinear Analysis of the C-Peptide Variable Related to Type 1-Diabetes Mellitus

Abstract

Type-1 diabetes mellitus is a chronic disease that is constantly monitored worldwide by researchers who are strongly determined to establish mathematical and experimental strategies that lead to a breakthrough toward an immunological treatment or a mathematical model that would update the UVA/Padova algorithm. In this work, we aim at a nonlinear mathematical analysis related to a fifth-order ordinary differential equations model that describes the asymmetric relation between C-peptides, pancreatic cells, and the immunological response. The latter is based on both the Localization of Compact Invariant Set (LCIS) appliance and Lyapunov’s stability theory to discuss the viability of implementing a possible treatment that stabilizes a specific set of cell populations. Our main result is to establish conditions for the existence of a localizing compact invariant domain that contains all the dynamics of diabetes mellitus. These conditions become essential for the localizing domain and stabilize the cell populations within desired levels, i.e., a state where a patient with diabetes could consider a healthy stage. Moreover, these domains demonstrate the cell populations’ asymmetric behavior since both the dynamics and the localizing domain of each cell population are defined into the positive orthant. Furthermore, closed-loop analysis is discussed by proposing two regulatory inputs opening the possibility of nonlinear control. Additionally, numerical simulations show that all trajectories converge inside the positive domain once given an initial condition. Finally, there is a discussion about the biological implications derived from the analytical results.

1. Introduction

Diabetes is a severe long-term condition ranking as one of the first ten causes of death in adults; according to global estimations, around four million people worldwide died in 2017 from this disease. Since the year 2000, the International Diabetes Federation (IDF) has reported the regional, national, and global occurrence of diabetes, indicating that the worldwide population with diabetes may increase from 463 to 700 million in the next two decades [1]. Diabetes mellitus (DM) is a long-term condition resulting from the inability of the pancreas to produce enough insulin (type 1 diabetes, T1D) or from the incapability of the pancreas to process the insulin that the body produces (type 2 diabetes, T2D) [2]. Thus, an increase in blood glucose leads to a non-symmetric behavior in the body over time.

Some complications related with high glucose blood include hypertension, kidney failure, lower limb amputation, nerve damage, stroke, and blindness [3]. Studies investigating the trends in diabetes prevalence have been conducted since 2000 [4], including the diabetes prevalence forecast for 2030 [5], 2035 [6], 2040 [4], 2045 [7], and 2060 [8], based on the national and regional data, where the results were overwhelming. Recently, several mathematical models have been published describing the process of glucose-insulin into the regulatory system, and the so-called Bergman’s Minimal Model is the most highlighted [9,10].

Recently, research applying mathematical algorithms to describe diabetes behavior and its outcomes is gaining attention [11,12]. Some of these algorithms are focused on both modeling insulin receptor or the body’s insulin-glucose dynamics, diabetes cost-effectiveness, and glucose tolerance testing [13]. However, models used to describe the glucose’s dynamic, insulin transport, and accuracy of glucose measurements, are challenging to assess in vivo.

Therefore, studying these non-symmetrical metabolic processes by mathematical approaches can help to understand these dynamics [14]. On the other hand, models based on Ordinary Differential Equations (ODEs) have been widely applied to describe real-life systems in physics, engineering, economics, and biomedicine. In particular, ODE models become a promising alternative to describe within-host dynamics, infectious or viral diseases, and even complex biomedical behaviors of the human body [15].

Currently, clinical studies and in silico data have demonstrated that C-peptide administration reduces renal disfunction, and combinations with insulin helps avoid microvascular issues. Hence, patients with C-peptide persistence are less prone to long-term complications than those without it [16]. The study of $\beta$ cell population dynamics in long-time intervals becomes a key to understand the prevalence of C-peptide secretion in T1D [17].

Some studies demonstrate that C-peptide levels drop exponentially in the first seven years after diagnosis and could continue dropping throughout the years at a slower rate. Nonetheless, log-transformed C-peptide levels permit establishing differences, both pathophysiological and immunological, between glucose and pancreatic cells, giving essential knowledge to understand $\beta$ cell survival. Therefore, broader attention should be paid to the progression of C-peptide loss in a longer duration of T1D, even with special focusing on the patient’s age [18].

The Localization of Compact Invariant Set (LCIS) method is a reliable method commonly used in nonlinear ODE models with mathematical-biological implications, see [19,20]. This method helps to provide sufficient or necessary conditions that lead to a broad understanding of the long-term behavior of a dynamical model, even to establish requirements for possible treatments or reduce some undesired cell populations proliferation [19,20,21].

In the particular case of T1D, the LCIS method permits analysis of the $\beta$ cell behavior in the presence of glucose [22] or with the immunological response [23]. The ODE’s mathematical model was initially presented in [24] describing the dynamics between cytotoxic T cells, the $\beta$ cell population, and the peptide level as a result of their interactions.

Our objective is to provide the conditions for a localizing domain, understand the global behavior of T1D’s cell dynamics, and give viable cells stabilization conditions. Hence, our hypothesis aims at how maximum population cells behave in time, based on upper bounds computations.

The organized sections of this work are presented as follows. The first section describes a general scheme for the fifth-order nonlinear mathematical model where upper bounds for all variable states hold when the positive orthant domain is satisfied by the nonlinear model’s positiveness. Some of the proposed localizing functions have no mathematical restriction on how they are defined or in the quantity limitations associated with a particular upper bound; however, the proposed function must not be the first integral, see [25,26]. Discussion resulting from applying local asymptotic stability by Lyapunov indirect method given the equilibrium point led to analyzing the stability criteria by closed-loop Lyapunov in which the input controls are analyzed. The second section shows some simulations that validate our previous mathematical results, and the last section presents the main conclusions of this research.

2. Preliminaries of Localization of Compact Invarian Sets Method

This section presents the necessary background to define the localizing domain that contains all the compact invariant sets of a nonlinear system represented by first-order ODEs. The general method of LCIS was proposed by Krishchenko and Starkov in [25,26] to determine the domain where all compact invariant sets of a differential equations system are located. This method is helpful to understand the long-time behavior of first-order ODE systems. This is considered an autonomous nonlinear system represented by:

$$\dot{x}=f\left(x\right),$$

(1)

where $x\in {\mathbb{R}}^n$, $f\left(x\right)={(f_1\left(x\right),\dots ,f_n\left(x\right))}^T$ is a differentiable vector field. Let $h\left(x\right)\in C^{\infty }\left({\mathbb{R}}^n\right)$ be a function such that h is not the first integral of the system (1). The function h is exploited in the solution of the localization problem of compact invariant sets and is called a localizing function. By ${h|}_U$ we denote the restriction of h on a set $U\subset {\mathbb{R}}^n$. $S\left(h\right)$ denotes the set $\{x\in {\mathbb{R}}^n\mid L_{f}h\left(x\right)=0\}$, where $L_{f}h\left(x\right)$ is the Lie derivative in the vector field of $f\left(x\right)$. In order to determine the localizing set, it is necessary to define $h_{inf}\left(U\right):=inf\{h\left(x\right)\mid x\in U\cap S\left(h\right)\}$ and $h_{sup}\left(U\right):=sup\{h\left(x\right)\mid x\in U\cap S\left(h\right)\}$.

Therefore, for any $h\left(x\right)\in C^{\infty }\left({\mathbb{R}}^n\right)$, all compact invariant sets of the system (1) located in U are contained in the set $K(U;h)$ defined as $\{x\in U\mid h_{inf}\left(U\right)\le h\left(x\right)\le h_{sup}\left(U\right)\}$, as well as, if $U\cap S\left(h\right)=\varnothing$, then there are no compact invariant sets located in U. Moreover, the Iterative Theorem can be applied to refine the localizing domain $K\left(h\right)$; this theorem is defined as follows [19,20,21,23]:

Theorem 1.

(Iterative Theorem). Let $h_m\left(x\right)$ be a sequence of $C^{\infty }$ differentiable functions where $m=0,1,2...$. Sets

$$K_0=K\left(h_0\right),K_m=K_{m-1}\cap K_{m-1,m},m>0,$$

with

$$\begin{array}{c}{K}_{m-1,m}=\left\{x\mid h_{m,inf}\le h_m\left(x\right)\le h_{m,sup}\right\},\hfill \\ h_{m,sup}=\underset{S\left(h_m\right)\cap K_{m-1}}{sup}{h}_m\left(x\right),\hfill \\ h_{m,inf}=\underset{S\left(h_m\right)\cap K_{m-1}}{inf}{h}_m\left(x\right),\hfill \end{array}$$

contain any compact invariant set of the system (1) and

$$K_0\supseteq K_1\supseteq \dots \supseteq K_m\supseteq \dots .$$

In summary, the general methodology to compute the LCIS of a nonlinear dynamical system described by first-order ODEs is as follows [27]:

• A localizing function denoted as $h\left(x\right)$ must be proposed. $h\left(x\right)$ is a function that can represent a specific shape, such as a plane, hyperplane, cylinder, or sphere; in terms of the system’s parameters and state variables.
• Computing the Lie derivative of $h\left(x\right)$, defined as $L_{f}h$.
• Calculate the infimum ($h_{inf}$) and supremum ($h_{sup}$) by computing $L_{f}h=0$. From the latter, two cases can result:

  (a)

  $S\left(h\right)$ is compact, the Lagrange multiplier method or the polytope approximation may be applied;

  (b)

  the sign of $S\left(h\right)$ cannot be defined, a mapping must be performed to determine the sign of $h\left(x\right){|}_{S\left(h\right)}$.

• If it is not possible to define the sign of $S\left(h\right)$, the localization problem is not yet solved; therefore, a new localizing function must be proposed, and the process is restarted.
• In the case of a satisfactory localizing domain, Theorem (1) could be applied to refine the localizing domain $K\left(h\right)$.

This methodology can be applied until a satisfactory solution is achieved.

3. Mathematical Model of Type-1 Diabetes Mellitus Related to C-Peptide

The mathematical model of Type-1 Diabetes Mellitus (T1DM) related to C-peptides was proposed by Mahay and Edelstein-Keshet, in 2007 [24], involving the immune response as the main factor that leads to a decrease in the $\beta$ cell population in the body. It consists of five first-order ODEs describing the dynamical interaction between activated T cells ($A\left(t\right)$), memory T cells ($M\left(t\right)$), effector T cells ($E\left(t\right)$), the C-peptide level ($p\left(t\right)$), and the $\beta$ cell ($B\left(t\right)$) populations at time t. Therefore, the T1DM model related to the C-peptide is given as follows:

$$\begin{array}{ccc}\hfill \frac{dA}{dt}& =& (\sigma +\alpha M)\frac{p^n}{k_1^n+p^n}-(\beta +{\delta }_A)A-ϵA^2,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \frac{dM}{dt}& =& \beta 2^{m_1}\frac{a{k}_2^m}{k_2^m+p^m}A-\frac{p^n}{k_1^n+p^n}\alpha M-{\delta }_{M}M,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \frac{dE}{dt}& =& \beta 2^{m_2}(1-\frac{a{k}_2^m}{k_2^m+p^m})A-{\delta }_{E}E,\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \frac{dp}{dt}& =& REB-{\delta }_{p}p,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill \frac{dB}{dt}& =& -\kappa EB;\hfill \end{array}$$

(6)

where Equations (2)–(4) correspond to the population level of activated, memory, and effector T cells; Equation (5) represents the peptide level and the remaining population of $\beta$ cells by Equation (6). The parametrization and units of the model’s equations are presented in

Table 1. Parameter description and units of T1DM related to C-peptides [24].

Parameter	Description	Value	Units
$\sigma$	Influence of naive T cells in the thymus	1–10	day${}^{-1}$
$\alpha$	Production rate of A per M	1–5	${\mathrm{day}}^{-1}$
$\beta$	Cell division rate	1–6	day${}^{-1}$
${\delta }_A$	Mortality index, activated T cells	≈0.01	day${}^{-1}$
${\delta }_M$	Mortality index, memory T cells	≈0.01	day${}^{-1}$
${\delta }_E$	Mortality index, effector T cells	$0.3$	day${}^{-1}$
${\delta }_p$	Peptide turnover rate	0–1	day${}^{-1}$
$ϵ$	Competition parameter T cell	1–5 $\times {10}^{-2}$	cells${}^{-1}$ day${}^{-1}$
$k_1$	Peptide level for $1/2$ of maximum
	activation	2	$peptide$$units$
$k_2$	Peptide level for $1/2$ of the maximum
	memory cells	1	$peptide$$units$
m	Hill’s coefficient, production of
	memory cells	2	−
n	Hill’s coefficient for activation of T cells	3	−
$2^{m_1}$	Maximum number of memory cells
	produced by T cells proliferating	8	−
$2^{m_2}$	Number of effector cells produced by
	proliferating T cells	60	−
a	Maximum fraction of memory cells
	produced	<1	−
R	Peptide accumulation rate	${10}^{-5}$	cells${}^{-1}$ day${}^{-1}$
$\kappa$	Death of $\beta$ cells by effector T cells	$0.14\times {10}^{-6}$	cells${}^{-1}$ day${}^{-1}$

.

Furthermore, to fulfill the positivity of the system (2)–(6), we evaluated each state variable at the border, i.e., $A=E=M=p=B=0$. Evaluating Equation (6), we obtained that $\frac{dB}{dt}=0$; Equation (5) gives that $\frac{dp}{dt}=REB$; Equation (4) gives that $\frac{dE}{dt}=\beta 2^{m_2}(1-\frac{a{k}_2^m}{k_2^m+p^m})A$; whereas Equation (3) gives that $\frac{dM}{dt}=\beta 2^{m_1}\frac{a{k}_2^m}{k_2^m+p^m}$; finally, from Equation (2), we obtained that $\frac{dA}{dt}=(\sigma +\alpha M)\frac{p^n}{k_1^n+p^n}$; allowing us to conclude that, given nonnegative initial conditions, the system’s dynamics are located in the non-negative orthant, i.e., they are located into the following domain:

$${\mathbf{R}}_{+,0}^5=\{A\left(t\right)\ge 0,M\left(t\right)\ge 0,E\left(t\right)\ge 0,p\left(t\right)\ge 0,B\left(t\right)\ge 0\}.$$

(7)

4. Localization of Compact Invariant Sets-Peptide Variable Analysis

Localizing the compact invariant domain of a dynamical system depends on the system’s complexity. Sometimes, it is possible to define the domain of attraction that contains all compact invariant sets by employing only one localizing function, resulting in symmetric shapes, such as ellipsoids, paraboloids, and cylinders [26]. These shapes are frequently obtained in three-dimensional systems.

However, biological systems are often modeled by more than three dimensions, making it impossible to define symmetric shapes by only one function and, at the same time, ensure all the dynamics of a system are bounded. Biological systems usually need more than one localizing function to describe the system’s variables’ maximum and minimum bounds [19]; therefore, the compact localizing domain is characterized by an asymmetric domain. Hence, the compact localizing domain and ultimate bounds for a T1DM related to C-peptide are achieved by exploring three localizing functions. First, we compute the maximum population of $\beta$ cells with the following localizing function

$$h_1=B-lnB,$$

(8)

whose Lie derivative is given by $L_f{h}_1=-\kappa EB-\left[\frac{-\kappa EB}{B}\right]$, defining the set $S\left(h_1\right)=\left\{L_f{h}_1=0\right\}$ as

$$S\left(h_1\right)=\left\{-\kappa EB+\kappa E=0\right\},$$

(9)

and, after solving for B, the set $S\left(h_1\right)$ is defined as

$$S\left(h_1\right)=\left\{B=1\right\};$$

(10)

further, expressing the constraint $h_1{|}_{S\left(h_1\right)}=B-lnB$, and substituting $S\left(h_1\right)$, the maximum value of the function $h_1$ is as follows

$$K\left(h_1\right)=\left\{h_1\le h_1{|}_{S\left(h_1\right)}:=1\right\}.$$

(11)

Therefore, the location set of the $\beta$ cell population is

$$K_B=\{B\left(t\right)\le B_{max}:=1\}.$$

(12)

Now, to determine the upper bound for the C-peptide level, the following localizing function is proposed

$$h_2=p-lnp+B,$$

(13)

where its Lie derivative is defined by $L_f{h}_2=REB-{\delta }_{p}p-\left[\frac{REB-{\delta }_{p}p}{p}\right]-\kappa EB$. The set is determined and analyzed by $S\left(h_2\right)=\left\{L_f{h}_2=0\right\}$, giving

$$S\left(h_2\right)=\left\{REB-{\delta }_{p}p-\frac{REB}{p}+{\delta }_p-\kappa EB=0\right\},$$

(14)

and it can be defined in terms of the interest variable of the localizing function as $S\left(h_2\right)=\left\{p=1-\frac{REB}{{\delta }_{p}p}-(\kappa -R)\frac{EB}{{\delta }_p}\right\}$, where, after some algebraic manipulation gives

$$S\left(h_2\right)=\left\{p=1\right\},$$

(15)

as long as the following condition can be satisfied

$$R\le \kappa .$$

(16)

It is important that the constraint can be expressed by $h_2{|}_{S\left(h_2\right)}=p-lnp+B$. Substituting (15) into $h_2{|}_{S\left(h_2\right)}$ and applying Theorem (1), the set $K\left(h_2\right)$ is defined as

$$K\left(h_2\right)\cap K\left(h_1\right)=\left\{h_2\le h_2{|}_{S\left(h_2\right)}:=1+B_{max}\right\};$$

(17)

allowing us to compute the set $K\left(h_2\right)$, which defines the maximum C-peptide level as

$$K_p=\{p\left(t\right)\le p_{max}:=2\}.$$

(18)

Finally, an upper bound for T cells is computed through the following localizing function

$$h_3=A+M+E,$$

(19)

whose Lie derivative is given by

$$\begin{array}{ccc}\hfill L_f{h}_3& =& \{\left(\sigma +\alpha M-\alpha M\right)\frac{p^n}{k_1^n+p^n}+(\beta 2^{m_2}-\beta -{\delta }_A-ϵA)A\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& & +\left(2^{m_1}-2^{m_2}\right)\frac{\beta a{k}_2^m}{k_2^m+p^m}A-{\delta }_{M}M-{\delta }_{E}E\};\hfill \end{array}$$

(21)

hence, after some algebraic rearrangement and mathematical analysis, the set $S\left(h_3\right)=\left\{L_f{h}_3=0\right\}$ is defined as

$$S\left(h_3\right)=\left\{A=\frac{\sigma p^n}{{\delta }_A(k_1^n+p^n)}+\frac{\left(\beta 2^{m_2}-\beta \right)}{{\delta }_A}A-\frac{ϵ}{{\delta }_A}{A}^2-\frac{{\delta }_M}{{\delta }_A}M-\frac{{\delta }_E}{{\delta }_A}E\right\},$$

(22)

as long as the next condition holds

$$2^{m_1}\le 2^{m_2}.$$

(23)

Now, substituting the previous results and some algebraic manipulation, the constrain $h_3{|}_{S\left(h_3\right)}=A+M+E$ is defined by

$$h_3{|}_{S\left(h_3\right)}=\frac{\sigma }{{\delta }_A}+\frac{{\left(\beta 2^{m_2}-\beta \right)}^2}{4ϵ{\delta }_A},$$

as long as the following conditions must be satisfied at all time

$$\begin{array}{ccc}\hfill {\delta }_A& \le & {\delta }_M,\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill {\delta }_A& \le & {\delta }_E;\hfill \end{array}$$

(25)

then, it is possible to define the set $K\left(h_3\right)$ as

$$K\left(h_3\right)=\left\{h_3\le h_3{|}_{S\left(h_3\right)}:=\frac{\sigma }{{\delta }_A}+\frac{{\left(\beta 2^{m_2}-\beta \right)}^2}{4ϵ{\delta }_A}\right\}.$$

Summarizing the results shown through this section, the following statement is formulated regarding the ultimate bounds for the T1DM related to the C-peptide system.

Theorem 2.

If the conditions (16), (24), (25) are fulfilled, all the compact invariant sets of the T1DM related to C-peptide system (2)–(6) lie within the following domain location

$$\begin{array}{c}\hfill Kse=K_B\cap K_p\cap K_A\cap K_M\cap K_E,\end{array}$$

(26)

where

$$\begin{array}{ccc}\hfill K_B& =& \{B\left(t\right)\le B_{max}:=1\};\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill K_p& =& \{p\left(t\right)\le p_{max}:=2\};\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill K_A& =& \left\{A\left(t\right)\le A_{max}:=\frac{\sigma }{{\delta }_A}+\frac{{\left(\beta 2^{m_2}-\beta \right)}^2}{4ϵ{\delta }_A}\right\};\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill K_M& =& \left\{M\left(t\right)\le M_{max}:=\frac{\sigma }{{\delta }_A}+\frac{{\left(\beta 2^{m_2}-\beta \right)}^2}{4ϵ{\delta }_A}\right\};\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill K_E& =& \left\{E\left(t\right)\le E_{max}:=\frac{\sigma }{{\delta }_A}+\frac{{\left(\beta 2^{m_2}-\beta \right)}^2}{4ϵ{\delta }_A}\right\}.\hfill \end{array}$$

(31)

The skewness correlation between C-peptides and cells was also demonstrated in [24], as the time scale of the peptide dynamics is faster (hours) than the time scale of the cell dynamics (days), and thus an almost steady state is assumed in the peptide. The C-peptide clinical test, which is widely applied to measure pancreatic $\beta$ cell function [28].

Considering the mathematical function $dp/dt=0$, leads to the variable C-peptide as $p=(REB/{\delta }_p)$. In this case, the state variable is far from being defined as an invariant plane in the mathematical scope; further, C-peptide represents a function that relays in the $\beta$ cell population with the immune response’s presence through effector cell populations. A disadvantage of analyzing the C-peptide in terms of other variables implies that the maximum carrying capacity of $\beta$ cells can be estimated in a general scheme. This research contributes by analyzing a whole model with the LCIS method to establish a scheme where the clinical test interpretation can lead to a mathematical preamble approach.

4.1. Local Stability

In this subsection is presented the mathematical results applying the Lyapunov indirect method and considering the equilibrium point $\left(A^{*},M^{*},E^{*},p^{*},B^{*}\right)=\left(0,0,0,0,0\right)$ in the positive orthant. To determine if the equilibrium is locally stable, the system of Equations (2)–(6) is linearized. First, the system’s Jacobian matrix (J) is defined as follows

$$J=\left[\begin{array}{ccccc}-(\beta +{\delta }_A)-2ϵA& \alpha \frac{p^n}{k_1^n+p_n}& 0& (\sigma +\alpha M){\varphi }_2& 0\\ \beta 2^{m_1}{\varphi }_3& -\frac{p^n}{k_1^n+p_n}\alpha -{\delta }_M& 0& {\varphi }_1{2}^{m_1}-{\varphi }_2\alpha M& 0\\ \beta 2^{m_2}(1-{\varphi }_3)& 0& -{\delta }_E& -{\varphi }_1{2}^{m_2}& 0\\ 0& 0& RB& -{\delta }_p& RE\\ 0& 0& -kB& 0& -kE\end{array}\right],$$

(32)

where

$$\begin{array}{ccc}\hfill {\varphi }_1& =& \beta \frac{am{k}_2^m{p}^{m-1}}{{\left(k_2^m+p^m\right)}^2}A,\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill {\varphi }_2& =& \frac{n{k}_1^n{p}^{n-1}}{{\left(k_1^n+p_n\right)}^2},\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill {\varphi }_3& =& \frac{a{k}_2^m}{k_2^m+p^m};\hfill \end{array}$$

(35)

evaluating matrix J at the equilibrium point, we obtained the expression

$$J\left(A^{*},M^{*},E^{*},p^{*},B^{*}\right)=\left[\begin{array}{ccccc}-(\beta +{\delta }_A)& 0& 0& 0& 0\\ \beta 2^{m_1}a& -{\delta }_M& 0& 0& 0\\ \beta 2^{m_2}(1-a)& 0& -{\delta }_E& 0& 0\\ 0& 0& RB& -{\delta }_p& 0\\ 0& 0& -kB& 0& 0\end{array}\right];$$

(36)

thus, the eigenvalues of (36) are ${\lambda }_1=-(\beta +{\delta }_A)$, ${\lambda }_2=-{\delta }_M$, ${\lambda }_3=-{\delta }_E$, ${\lambda }_4=-{\delta }_p$, and ${\lambda }_5=0$.

Since ${\lambda }_5=0$, it is impossible to conclude local stability for the equilibrium by applying the Lyapunov indirect method theorem in Equation (36). In summary, the system of Equations (2)–(6) has only one equilibrium point; therefore, local asymptotic stability is not evident. Hence, the design criteria in which the authors initially based the system (2)–(6) in [24] opens the possibility of considering control inputs to define a complementary model.

However, implementing the LCIS method provided a positive domain where all nonlinear system’s trajectories were held without a linear scheme or numerical approach; thus, establishing a solution to the system by defining the upper bounds given the conditions (16), (24), and (25). The domain defined by (26) contains the cell population; however, it is considered asymmetric regarding each cell dynamic.

Closed-Loop Analysis via Lyapunov Stability Criteria

In the particular case of biological models, proposing control inputs are complex to determine unless a real-world known variable can be measured or supplied in a laboratory, such as insulin. In this work, insulin is not directly involved; instead, we assumed that a more comprehensive understanding of blocking a direct targeting of the effector cells to the pancreatic cells would lead to unnecessary antigen scheme behavior.

Recent research suggests that a more in-depth development of the insulin proliferation due to the $\beta$ cell behavior. In [29], the authors concluded that researchers worldwide must continue monitoring T1D incidence trends. In contrast, research associated with prevention areas, early detection, and improved TID treatment continues. Furthermore, in [30], the authors tackled the use of protein biomarkers associated with risk factors in developing cardiovascular diseases when diabetes family antecedents prevail and pass in offspring from the gestational diabetes stage. They concluded that a deeper understanding of a leading cause that diabetes develops could improve this research topic.

Therefore, considering the system dynamic and the obtained previous results, we decided to analyze the system in a closed-loop scheme, proposing control inputs that guarantee its overall stability. The model described by Equations (2)–(6) is expressed as follows:

$$\begin{array}{ccc}\hfill \frac{dA}{dt}& =& (\sigma +\alpha M)\frac{p^n}{k_1^n+p^n}-(\beta +{\delta }_A)A-ϵA^2+U_1,\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill \frac{dM}{dt}& =& \beta 2^{m_1}\frac{a{k}_2^m}{k_2^m+p^m}A-\frac{p^n}{k_1^n+p^n}\alpha M-{\delta }_{M}M,\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill \frac{dE}{dt}& =& \beta 2^{m_2}(1-\frac{a{k}_2^m}{k_2^m+p^m})A-{\delta }_{E}E+U_2,\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill \frac{dp}{dt}& =& REB-{\delta }_{p}p,\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill \frac{dB}{dt}& =& -\kappa EB;\hfill \end{array}$$

(41)

where $U_1$ and $U_2$ represent the control inputs that could regulate the T cell population increase rate. To determine the criteria for each input, we propose the following candidate Lyapunov function

$$\begin{array}{c}\hfill V=q_{1}A+q_{2}M+q_{3}E+q_{4}p+q_{5}B,\end{array}$$

(42)

where its derivative is given by $\dot{V}=q_1\dot{A}+q_2\dot{M}+q_3\dot{E}+q_4\dot{p}+q_5\dot{B}$, with $q_1$, $q_2$, $q_3$, $q_4$, and $q_5$ as free parameters, after substituting Equations (37)–(41) into the derivative gives

$$\begin{array}{ccc}\hfill \dot{V}& =& q_1\left[\sigma \frac{p^n}{k_1^n+p_n}+\alpha M\frac{p^n}{k_1^n+p_n}-(\beta +{\delta }_A)A-ϵA^2+U_1\right]\hfill \\ & & +q_2\left[\beta 2^{m_1}\frac{a{k}_2^m}{k_2^m+p^m}A-\frac{p^n}{k_1^n+p_n}\alpha M-{\delta }_{M}M\right]\hfill \\ & & +q_3\left[\beta 2^{m_2}A-\frac{a{k}_2^m}{k_2^m+p^m}\beta 2^{m_2}A-{\delta }_{E}E+U_2\right]+q_4\left[REB-{\delta }_{p}p\right]+q_5\left[-\kappa EB\right].\hfill \end{array}$$

(43)

Now, analyzing the Equation (43), we concluded that $U_1$ and $U_2$ are able to establish stability conditions; therefore, $U_1$ and $U_2$ are defined as

$$\begin{array}{ccc}\hfill U_1& =& -\sigma \frac{p^n}{k_1^n+p_n}-\alpha M\frac{p^n}{k_1^n+p_n},\hfill \end{array}$$

(44)

$$\begin{array}{ccc}\hfill U_2& =& -\beta 2^{m_2}A,\hfill \end{array}$$

(45)

as long as the condition (46) holds

$$2^{m_1}<2^{m_2},$$

(46)

with

$$\begin{array}{ccc}\hfill q_1& =& q_2=q_3=q_4=1,\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill q_5& <& \frac{R}{\kappa },\hfill \end{array}$$

(48)

and, in order to guarantee asymptotically stability by the Lyapunov direct method, also the following inequality must hold

$$ϵ\left(A+\frac{(\beta +{\delta }_A)}{2ϵ}\right)+{\delta }_{M}M+{\delta }_{E}E+{\delta }_{p}p>\frac{{(\beta +{\delta }_A)}^2}{4ϵ}.$$

(49)

In summary, the condition for $q_5$, inequality (48), implies that, given the Equation (16), when $R=\kappa$, then $q_5<1$, in comparison with the positive free parameters (47) that are equal to one. Meanwhile, condition (46) holds as condition (23). This implies that set $K\left(h_3\right)$ in Equation (26) encompasses the system (37)–(41) only when $R=\kappa$; thus, Equation (49) is also contained in the positive domain of $K\left(h_3\right)$; leading us to a mathematical preamble that the system (2)–(6) is a baseline model that can guide a mathematical revaluation, i.e., a model where cell populations could be modified, in view of a possible treatment.

5. Numerical Simulations

This section presents numerical simulations obtained with the LCIS method. Figure 1, shows the behavior of the activated, effector, and memory T cells, as well as the behavior of the population level of $\beta$ cells and C-peptides. The parameters considered were those corresponding to Table 1, and the initial conditions were $A\left(0\right)=0.5$, $M\left(0\right)=0$, $E\left(0\right)=1$, $p\left(0\right)=0$, and $B\left(0\right)=1$. Figure 1 shows the number of circulating cells (scaled) against time (days); $A\left(t\right)$ is expressed as $\times {10}^3$ cells. $M\left(t\right)$($\times {10}^4$), $E\left(t\right)$($\times {10}^6$), $p\left(t\right)$ tends to be a small population of cells, and $B\left(t\right)$ is a fraction of the remaining $\beta$ cell mass. The $\beta$ cell population decreased by 40% during the immune attack. Since the model does not address the replenishment of the $\beta$ cells by reproduction or stem cell differentiation, the $\beta$ cell mass remains constant after this isolated immune response, [31]. The proposed initial conditions leading to the immune response was resolved without chronic disease or cyclic waves.

In Figure 2, we present a first approach of the upper bound for the variable $A\left(t\right)$, only if the conditions (24) and (25) are fulfilled; whereas, the immunological response in the presence of $\beta$ cell behavior is presented in Figure 3. Effector and memory cell dynamics are under the upper bound set $K\left(h_3\right)$, implying that C-peptide has a direct impact on them; therefore, in mathematical sense, the model has a proximity approach to the biological scheme. Clinical procedures need to be considered to ensure a reliable approach between the mathematical and the physical dynamics.

Figure 3 and Figure 4 show the maximum value of the set $K\left(h_3\right)$ for the two types of T cells. The condition of ${\delta }_A$, in (25), implies that the death rate of the memory T cells must be lower than the death rate by effector cells. Figure 5 presents the dynamics under the upper bound $K_p$, given by the localizing function $h_2$ and satisfying the condition (16). The secretion of the peptide is directly associated with the activation of T cells. When the C-peptide reaches high levels, memory cell production stops, and, consequently, the C-peptide is gradually cleared. High T cell levels are associated with an immune response to attack the $\beta$ cell population, while the C-peptide attempts to avoid their destruction.

Using the LCIS method, we determined the maximum $\beta$ cell population; therefore, when the $\beta$ cells are at the maximum, then the C-peptide level is at the minimum as long as the T cells remain inactivated, see Figure 1. However, the C-peptide secretion stops when the $\beta$ cells are gone; otherwise, its secretion remains active and waiting for the following $\beta$ cell-level change. An increased incidence of microvascular complications are correlated with low C-peptide levels. It would be interesting to determine whether C-peptide concentrations are associated with increased macrovascular morbidity and mortality. Moreover, the maximum population of $\beta$ cells is given by the set $K_B$, see Figure 6.

6. Conclusions

The localizing compact invariant set method provides the mathematical preamble to define the bounded positive invariant domain, i.e., the domain where all trajectories of the cell populations involved in T1DM are contained. It was also possible to mathematically describe the C-peptide level by proposing linear type localizing functions. The particularity in which the mathematical model is presented in [24], and discussed in this research implies that it represents a feasible scheme to analyze $\beta$ cell targeting by the immune response. Thus, the estimated numerical values in Table 1 hold a reliable approach that can lead to a deeper pancreatic cell population understanding for experimental research in the future.

C-peptides are a useful indicator of $\beta$ cell function, allowing discrimination between insulin-sufficient and insulin-deficient individuals with diabetes. Potential future uses of C-peptide are broad, including aiding appropriate diagnosis, guiding therapy choices, and predicting morbidity in diabetes; hence, the set of Equations (2)–(4) is one of the first nonlinear models involving a variable for C-peptide, and our results aim to contribute to future research involving a mathematical preamble.

The local stability of the systems through linearization was not concluded, since, in both cases, a matrix with a null eigenvalue was obtained, that is, one of the eigenvalues is equal to zero. Thus, this indicates the possibility of needing control inputs to ensure that the system regulates and breaks even.

The mathematical analysis of closed-loop systems suggests two control inputs directly related to the population of activated T cells and effector T cells. The control input $U_1$, see condition (44), implies the existence of a counterpart that prevents an increase in the population of activated T cells under the presence of $\beta$ cells by suppressing the C-peptide level and the number of memory T cells produced by the body. On the other hand, the control input $U_2$, see (45), is associated with the effector T cell population’s level, suggesting a population reduction effect of activated T cells to proliferate.

Therefore, a mathematical analysis considering control inputs based on a closed-loop system provides a theoretical basis to implement an immunotherapy treatment, if and only if, the conditions (46), (47), and (48) hold and, as a consequence, the condition (49) is also satisfied. In other words, these control inputs permit the conduction of all the cell populations to the desired state of equilibrium, i.e., being in symmetry with the desired level of each cell population.

This work did not discuss the idea of a nonlinear controller design at the moment; however, this is considered as future work given the conditions (44) and (45). We also intend to carry out the design of observers. We assumed, that the mathematical purpose of the observer is to identify or estimate those model’s variables for feedback and to implement it as a possible or feasible treatment. Since the model deals with cell populations that do not have an easy way to measure themselves, considering their development outside the body is still a goal for the future.