Dynamical Models in Intrinsic Lyapunov Methods

Ifeoma D. Omoko^1*, Oluwabusayo Oni², Taiwo G. Alabi³, Muritala B. Taoreed⁴

¹Department of Mathematics, Lagos State University, Ojo, Lagos, Nigeria.
²Department of Mathematics, Oregon State University, Oregon, USA.
³Department of Mathematics, University of Tennessee, USA.
⁴Department of Mathematics, Boise State University, Boise, Idaho, USA.

DOI: 10.4236/jamp.2025.134073                

Abstract

This study examines the construction of Intrinsic Lyapunov functions to investigate the asymptotic stability of equilibrium states analyzed under specific conditions involving interaction coefficients and compartmental transitions. It involves the qualitative properties of stratified compartmental models for drug user populations within heterogeneous communities using intrinsic Lyapunov methods. These models categorize populations into different compartments (e.g., susceptible individuals, drug users, and those undergoing treatment) and account for interactions among these subpopulations. The system is governed by nonlinear ordinary differential equations, with parameters representing initiation, recovery, and relapse rates. Additionally, the boundedness of solutions is established, ensuring that population variables (e.g. the number of drug users) remain within biologically feasible limits over a period of time. These findings provide theoretical insights into the long-term dynamics of drug user populations, highlighting the impact of prevention, treatment, and socio-environmental factors. The results underscore the effectiveness of intrinsic Lyapunov methods in understanding stratified population dynamics and guiding targeted intervention strategies to control drug prevalence.

Keywords

Nonlinear Differential Equations, Stability, Boundedness, Intrinsic Lyapunov Method, Population Models

1. Introduction

This research focuses on developing Lyapunov functions for systems of nonlinear differential equation by intrinsic Lyapunov methods to establish stability and boundedness in drug user population models. Unlike conventional Lyapunov functions, intrinsic Lyapunov functions are coordinate-independent and often tied to the system’s geometry or physical principles [1]. These functions are particularly useful in systems governed by conservation laws, Hamiltonian structures, or other intrinsic dynamics [2]. As dynamical systems theory has advanced, researchers have sought Lyapunov functions expressed in intrinsic terms, often utilizing differential geometry and functional analysis to describe stability properties independent of specific coordinate systems [3]. Intrinsic Lyapunov functions find applications across diverse fields such as in Control Systems—Stabilizing controllers for robotic systems and power networks, Mechanical Systems, Thermodynamic and Biological Systems. Research on intrinsic Lyapunov Methods are considerably few. Chin [4] treated the construction of the generalized intrinsic Lyapunov function and its extension for $\left(k+1\right)th,kth$ equations of a set of n first order nonlinear differential equations:

${\dot{x}}_{k+1}=f_{k+1}\left(x_1,x_2,x_3,\cdots ,x_n\right)$

${\dot{x}}_k=f_k\left(x_1,x_2,x_3,\cdots ,x_n\right)$

Shi-Zhang [5] examined qualitative properties such as stability and boundedness solutions of the non linear systems:

$x^{iv}+f_1\left(\ddot{x}\right)\stackrel{⃛}{x}+f_2\left(\dot{x},\ddot{x}\right)+f_3\left(\dot{x}\right)+f_4\left(\dot{x},\ddot{x}\right)=P\left(t,x,\dot{x},\ddot{x},\stackrel{⃛}{x}\right)$

where $f_1,f_2,f_3,f_4$ are continuous functions. Other approaches to intrinsic Lyapunov methods have been explored by Chin, Tunc [6]-[8]. In this research, the construction of the Lyapunov functions is carried out using population models of drug users. The dynamics of drug user populations in a heterogeneous population can be analyzed using mathematical models that incorporate stability theory. These models often describe how the number of drug users changes over time, influenced by intra- and inter-subpopulation interactions, as examined by Fikiri and Andreas [9]-[11]. Resmawan [12] examined the dynamic model analysis of the spread of drug addicts with educational effects on the model

$\frac{\text{d}N}{\text{d}t}=\text{Λ}-\mu N-{\text{Γ}}_{1}U-{\text{Γ}}_2{U}_e$

where $N=U+S+R+S_e+U_e+R_e$ is the total number of population and the variables are defined in their respective arguments. Other results on dynamical analysis include results from Oni, Olutimo, Omoko ([14]-[15] etc.).

2. Statement of the Problem

Our major emphasis in this research is the construction of the Lyapunov function using intrinsic methods. This method has advantages over explicitly methods since it overcomes the barriers and complexities involved in dealing with nonlinear systems of several compartmental variables. Lack of straight forward procedure in the Lyapunov function construction is minimized. The application of this method is a fundamental tool for analyzing the stability and boundedness of solutions to nonlinear system describing the dynamics of stratified population models.

3. Main Results

Mathematical Models

The population models are categorized into four classes: those at risk of using drugs(susceptible) denoted by S, lighter users L, Heavy users H and drug users on treatment T. Building on results of Shi-Zhang, Chin, Chin, Tunc, Fikiri ([4]-[9], Resmawan, [12]), the stability and boundedness of solutions for a stratified drug user model are examined through the following system:

$\begin{array}{l}\dot{S}=\delta -aS-bS+c_{1}L\\ \dot{L}=as+c_{2}H-d_{1}L\\ \dot{H}=eL+c_{3}T-d_{2}H\\ \dot{T}=fH-d_{3}T\end{array}$ (1)

Each compartment represents a category within the population. S, Susceptible individuals are at risk of using drugs. L which are light users have started using drugs but are not heavily dependent. H, are heavy users with severe drug dependence and consequently, T, individuals in treatment, are attempting to recover. Parameters are; $\delta$ —total population, b—natural death rate to the susceptible class; a—drug mortality rate of the susceptible class; $c_1$ —light users rate at first exposure to the susceptible class, $d_1$ drug mortality rate among light users, $c_2$ —heavy users first exposure to the susceptible class, e—migration rate from heavy users to light users due to effective treatment, $c_3$ quality treatment to improve recovery rates, $d_2$ drug mortality rate on the heavy users, f—rate of transition from light users to heavy users due to treatment relapses and $d_3$ , Relapse rate among individuals in treatment.

4. Construction of Intrinsic Lyapunov Function

This method is for a general class of non-linear systems expressed in state variables as n first order non-linear equations. It applies integration by parts procedure, derives a Lyapunov function directly from the differential equation under study. For this, the integration is along trajectories and the limits for the integral with respect to time are from zero to t, [4]. The derivation of the Lyapunov function V and its derivative, $\dot{V}$ are based on the equation

$V+{\displaystyle {\int }_0^t-\dot{V}\text{d}\tau }=0$

and do not require the gradient of the scalar function V to be obtained [6] [16].

Definition 1. An intrinsic Lyapunov function $V\left(x\right)$ is a scalar function defined globally (or invariantly with respect to the system’s structure), satisfying:

• Positive Definiteness: $V\left(x\right)>0$ for all $x$ , and $V\left(x\right)=0$ .

• Monotonic Decrease: Along trajectories of the system, the derivative $\dot{V}\left(x\right)=\frac{\text{d}}{\text{d}t}V\left(x\right)\le 0$. This ensures $V\left(x\right)$ decreases or remains constant over time.

From the systems (1) we derive the following differentials as

$\frac{\text{d}S}{\text{d}L}=\frac{\delta -aS-bS+c_{1}L}{as+c_{2}H-d_{1}L}$

$\frac{\text{d}L}{\text{d}H}=\frac{as+c_{2}H-d_{1}L}{eL+c_{3}T-d_{2}H}$

$\frac{\text{d}H}{\text{d}T}=\frac{eL+c_{3}T-d_{2}H}{fH-d_{3}T},$

$\frac{\text{d}H}{\text{d}T},\frac{\text{d}S}{\text{d}T},\frac{\text{d}S}{\text{d}H},\frac{\text{d}L}{\text{d}T}$ , etc., and obtain the integrals by parts as follows;

$\begin{array}{c}{I}_1=\frac{a{S}^2}{2}+c_{2}HS+\delta L-\left(a+b-d_1\right)LS+\frac{c_1{L}^2}{2}\\ -{\displaystyle {\int }_0^t[-a{d}_1{S}^2-\left(a{c}_1+b{c}_1\right)L^2-\left(a^2+2ab+b^2-d_1^2+c_{2}e\right)LS}\\ -\left(c_2{d}_2-c_2{d}_1\right)HS-c_2{c}_{3}ST-\left(a+b\right)\delta L]\text{d}t\end{array}$

$\begin{array}{c}{I}_2=e{L}^2+c_{3}LT-\left(d_2+d_1\right)HL-aHS-c_2{H}^2\\ -{\displaystyle {\int }_0^t[\left(c_{3}f+d_1^2+d_2^2-a{c}_1\right)HL}-\left(c_3{d}_3+c_2{d}_3\right)LT\\ -d_{2}e{L}^2+\left(a^2+ab+b^2-a{d}_1\right)HS-c_2{d}_1{H}^2-a\delta H]\text{d}t\end{array}$

$\begin{array}{c}{I}_3=f{H}^2-\left(d_3+d_2\right)HT-eLT-c_3{T}^2-{\displaystyle {\int }_0^t[-d_{3}f{H}^2+c_3{d}_3{T}^2}\\ -\left(d_3^2+c_{2}e+d_2{d}_3\right)HT+\left(e{d}_3+e{d}_2\right)LT-aeST]\text{d}t\end{array}$

$\begin{array}{c}{I}_4=eLS-c_{3}ST-d_{2}HS-\delta H+\left(a+b\right)HS+c_{1}LH\\ -{\displaystyle {\int }_0^t[\left(ae+a{c}_1\right)S^2}+\left(c_{2}e-c_{3}f-d_2^2-a^2-2ab-b^2+c_1{c}_2\right)HS\\ -\left(e{d}_1+e{d}_2+c_1{d}_1\right)LS+\left(c_3{d}_3-c_3{d}_3\right)ST-\left(a{c}_1+b{c}_1\right)HL]\text{d}t\end{array}$

$\begin{array}{c}{I}_5=fHS-\left(a+b-d_3\right)ST-\delta T+c_{1}LT-{\displaystyle {\int }_0^t[efLS+\left(c_{3}f+d_3^2-a^2-ab\right)ST}\\ -\left(d_{2}f+d_{3}f\right)HS+\left(a+b\right)\delta T+\left(a{c}_1+b{c}_1-c_1{d}_1\right)LT]\text{d}t\end{array}$

$\begin{array}{c}{I}_6=fHL-\left(d_3-d_1\right)LT-aST-c_{2}HT-{\displaystyle {\int }_0^t[ef{L}^2+\left(c_{3}f+d_3^2-a{c}_1-c_{2}e-d_1^2\right)LT}\\ -\left(d_{2}f+d_{3}f\right)HL-a\delta T+\left(a^2+ab+a{d}_1\right)ST-c_2{c}_3{T}^2+\left(c_2{d}_2+c_2{d}_1\right)HT]\text{d}t\end{array}$

From $I_i,i=1,\cdots ,6$ , the Lyapunov function are the terms outside the integrals given as

$\begin{array}{c}2V=a{S}^2+2\left(c_2-a-d_{2}a+b\right)HS+\delta \left(L-H+T\right)+\left(d_1-a-b+e\right)LS\\ +\left(c_1+e\right)L^2+\left(c_3-e+c_1-\left(d_3-d_1\right)\right)LT+(f+c_1-\left(d_2+d_1\right)HL\\ +\left(f-c_2\right)H^2-(c_2-\left(d_3+d_2\right)HT-\left(a+b-d_3+c_3+a\right)ST-c_3{T}^2\end{array}$ (2)

Time derivatives along system trajectories (1) are the values within the integrals which do appear as

$\begin{array}{c}\dot{V}=-a{d}_1{S}^2-\left(a{c}_1+b{c}_1\right)L^2-\left(a^2+2ab+b^2-d_1^2+c_{2}e\right)LS\\ -\left(c_2{d}_2-c_2{d}_1\right)HS-c_2{c}_{3}ST-\left(a+b\right)\delta L+\left(c_{3}f+d_1^2+d_2^2-a{c}_1\right)HL\\ -\left(c_3{d}_3+c_2{d}_3\right)LT-d_{2}e{L}^2+\left(a^2+ab+b^2-a{d}_1\right)HS-c_2{d}_1{H}^2-a\delta H\\ -d_{3}f{H}^2+c_3{d}_3{T}^2-\left(d_3^2+c_{2}e+d_2{d}_3\right)HT+\left(e{d}_3+e{d}_2\right)LT-aeST\\ +\left(ae+a{c}_1\right)S^2+\left(c_{2}e-c_{3}f-d_2^2-a^2-ab-ab-b^2+c_1{c}_2\right)HS\\ -\left(e{d}_1+e{d}_2+c_1{d}_1\right)LS+\left(c_3{d}_3-c_3{d}_3\right)ST-\left(a{c}_1+b{c}_1\right)HL+efLS\\ +\left(c_{3}f+d_3^2-a^2-ab\right)ST-\left(d_{2}f+d_{3}f\right)HS+\left(a+b\right)\delta T\\ +\left(a{c}_1+b{c}_1-c_1{d}_1\right)LT+ef{L}^2+\left(c_{3}f+d_3^2-a{c}_1-c_{2}e-d_1^2\right)LT\\ -\left(d_{2}f+d_{3}f\right)HL-a\delta T+\left(a^2+ab+a{d}_1\right)ST-c_2{c}_3{T}^2+\left(c_2{d}_2+c_2{d}_1\right)HT\end{array}$

This method provides a geometric perspective on stability theory and boundedness analysis describing the dynamics of drug users, rehabilitated individuals, etc. It proves whether the state is asymptotically stable or the conditions which drugs prevalence stabilize or grow.

5. Stability Analysis

In the context of the Intrinsic Lyapunov Methods, stability determines whether small perturbations or changes in initial conditions lead to bounded or unbounded deviations in the solutions over time [17]. The stability of solutions to the nonlinear four-compartment drug model can be rigorously analyzed with the choice of the Lyapunov function, $V\left(S,L,H,T\right)$ to equilibrium, its time derivative $\dot{V}\left(S,L,H,T\right)$ and evaluating its behavior along the trajectories of the system [18]. This equilibrium represents the steady-state population in each compartment. A simple definition of stability is explained bellow.

Definition 2. The zero solution of Equation (1) is STABLE if given $ϵ>0$ , $t_0>0$ , there exists a $\nu \left(t_0,ϵ\right)>0$ such that whenever,

$\left|X_0\right|<\nu \left(t_0,ϵ\right)$ , $\left|X\left(t,t_0\right)\right|<ϵ$ for all $t\ge t_0$ .

Drawing from this definition in Omoko [15], we present the assumptions that stability theory are be proved.

Theorem 1. In addition to the assumptions imposed on the model parameters, $a,b,c_2,c_3,d_1,d_2,d_3,e,f$ in (1) as non negative values, this study also assumed that

$\left(a\left(a+b\right)+a\left(a+b+e-d_1\right)\right)>0,\left(\left(a+b\right)d_1{d}_2{d}_3\right)>0,\left(c_3-c_2-d_1-d_3\right)>0$

$\left(c_2-c_3-d_2-d_3\right)>0,\left(a+b-d_1\right)>0$ ,

then the zero solutions to (1) satisfy $\left|S\right|\to 0$ , $\left|L\right|\to 0$ , $\left|H\right|\to 0$ , $\left|T\right|\to 0$ as $t\to \infty$ .

Proof: We select (2) as the Lyapunov function which can be rewritten as

$\begin{array}{c}2V=k_1{S}^2+k_2{L}^2+k_3{H}^2+k_4{T}^2+k_{5}LS+k_6{\left(H+L\right)}^2\\ +k_7{\left(L+T\right)}^2+k_8{\left(H+S\right)}^2+k_9{\left(S+T\right)}^2+k_{10}{\left(H+T\right)}^2\end{array}$ (3)

$k_1=\left(4a+b-c_2+c_3-d_1+d_2-d_3-e\right)>0$

$k_2=\left(-a-b-3c_1-c_3+d_1+d_2+d_3+e-f\right)>0$

$k_3=\left(-b-c_1-c_2+d_1+d_2-d_3\right)>0$

$k_4=\left(2a+b-c_1-c_2-c_3-\delta -d_1+d_2+d_3+e\right)>0$

$k_5=\left(d_1-a-b+e\right)>0,k_6=\left(f+c_1-\left(d_2+d_1\right)\right)>0$

$k_7=\left(c_3-e+c_1-\left(d_3-d_1\right)\right)>0,k_8=\left(c_2-a-d_{2}a+b\right)>0$

$k_9=\left(a+b-d_3+c_3+a\right)>0,k_{10}=\left(f-c_2\right)H^2-c_2-\left(d_3+d_2\right)>0$

Therefore, $V\ge 0$ along system trajectories (1). In summary,

$\begin{array}{c}\dot{V}=-[{\tau }_1{S}^2+{\tau }_2{L}^2+{\tau }_3{H}^2+{\tau }_4{T}^2+{\tau }_{5}LS\\ +{\tau }_{6}HL+{\tau }_{7}LT+{\tau }_{8}HS+{\tau }_{9}ST+{\tau }_{10}HT]\end{array}$

where

${\tau }_1=a\left(a+b\right)+a\left(a+b+e-d_1\right),$

${\tau }_2=\left(a+b-e-d_1\right)c_1+d_1\left(e-c_1\right)+e\left(d_1+d_2-c_1-f\right),$

${\tau }_3=d_2\left(c_2-f\right),{\tau }_4=c_3\left(d_2-c_2\right)+d_3\left(c_3-c_2-d_1-d_3\right)$

${\tau }_5={\left(a+b\right)}^2+e\left(a+b\right)+c_1\left(a+b\right)+a\left(e+c_1\right)+d_1\left(b+e\right)+e\left(d_2-c_2\right)$

${\tau }_6=c_1\left(b+c_2-d_3\right)+c_2\left(e+c_1\right)+d_2\left(c_1-d_1-d_2\right)+f\left(e-c_1-c_3-d_1+d_3\right)$

$\begin{array}{c}{\tau }_7=\left(a+b+d_3\right)c_1+d_1\left(e-c_1-c_3+d_3\right)+c_3\left(f+c_1-d_1-d_2\right)\\ +e\left(c_2-d_2-d_3+d_3\left(e-c_1-c_3-d_1+d_3\right)\right)\end{array}$

$\begin{array}{c}{\tau }_8=c_2\left(a+b\right)+d_2\left(a+b\right)+c_2\left(a+b-d_1\right)+d_2\left(b-c_2+d_2\right)\\ +f\left(2a+b+c_3+d_3\right)\end{array}$

$\begin{array}{c}{\tau }_9={\left(a+b\right)}^2+a\left(a+b\right)+c_3\left(a+b\right)+a\left(c_1+c_3+d_1-d_3\right)\\ +c_3\left(b-d_2\right)+d_3\left(2a+b+c_3-d_3\right)\end{array}$

${\tau }_{10}=c_2\left(e-c_3-d_1+d_3\right)+d_2\left(c_2-d_2-d_3\right)+f\left(c_2-c_3-d_2-d_3\right)$

$\dot{V}\le 0$

provided ${\tau }_i,i=1,\cdots ,10>0$ . Since $V>0,\dot{V}\le 0$ , the conditions ensure asymptotic stability. A well-posed model requires that the total population remains within realistic bounds.

6. Boundedness of Solutions

Boundedness of solutions for the nonlinear system describing the dynamics of a population ensures that the population within each stratified class does not grow indefinitely and remains biologically feasible and demonstrate how targeted interventions can stabilize or reduce heavy drug users. Introducing $P\left(t,S,L,H,T\right)$ into the system (1), the population models becomes

$\begin{array}{l}\dot{S}=\delta -aS-bS+c_{1}L\\ \dot{L}=as+c_{2}H-d_{1}L\\ \dot{H}=eL+c_{3}T-d_{2}H\\ \dot{T}=fH-d_{3}T+P\left(t,S,L,H,T\right)\end{array}$ (4)

$P\left(t,S,L,H,T\right)$ is a functional response [19], [20], which involves all the four compartments. It include creating total awareness to the susceptible class, drug users, of dangers and negative effect on hard drugs, highlighting effective invention strategies. This is suitable for creating boundedness of the Lyapunov function $V\left(S,L,H,T\right)$ to ensure the drug population do not grow endemically and indefinitely. From the definition on boundedness.

Definition 3. A solution $X\left(t,S,L,H,T\right)$ of Equation (4) is said to bounded if exists a constant $K>0$ such that $\left|X\left(t,S,L,H,T\right)\right|<K$ , where $K$ may depend on each solution.

Theorem 2. In addition to the assumptions imposed on the parameters in Theorem 1, let $\left|P\left(t,S,L,H,T\right)\right|\le \mu \left(t\right)\left(\left|S\right|+\left|L\right|+\left|H\right|+\left|T\right|\right)$ where $\mu$ is a positive continuous function of t, ${\displaystyle {\int }_0^t\mu \left(s\right)\text{d}s}\le \Omega <\infty$ , $\Omega >0$ , a constant, then there exist another constant, $\Gamma >0$ such that every solution of (4) satisfies $\left|S\right|<\Gamma$ , $\left|L\right|<\Gamma$ , $\left|H\right|<\Gamma$ , $\left|T\right|<\Gamma$ , $\forall t>0$ .

Proof. Define an intrinsic Lyapunov function $V\left(S,L,H,T\right)$ (2) as the total population’s deviation from a biologically realistic bound. Automatically the Lyapunov function is positive definite from (3). Time derivative with $P\left(t,S,L,H,T\right)$ using notation

$\dot{V}=\frac{\partial V}{\partial S}\dot{S}+\frac{\partial V}{\partial L}\dot{L}+\frac{\partial V}{\partial H}\dot{H}+\frac{\partial V}{\partial T}\dot{T},$

where

$\begin{array}{c}\frac{\partial V}{\partial S}\dot{S}=(h\left(a-a+b+c_2-d_2\right)-t\left(a+a+b+c_3-d_3\right)\\ +l\left(-a-b+d_1+e\right)+as)\left(\delta -aS-bD+c_{1}L\right)\\ =-[a\left(a+b\right)S^2+\left({\left(a+b\right)}^2+e\left(a+b\right)+c_1\left(a+b\right)\right)SL\\ +\left({\left(a+b\right)}^2+a\left(a+b\right)+c_3\left(a+b\right)\right)ST+\delta (bH+aS\left(a+b-e\right)L\\ +\left(2a+b\right)T+\left(\left(a+b+d_3\right)c_1\right)LT+\left(\left(a+b-e-d_1\right)c_1\right)L^2\\ +\left(c_1\left(b+c_2-d_3\right)\right)HL+\left(c_2\left(a+b\right)+d_2\left(a+b\right)\right)HS\end{array}$

$\begin{array}{c}\frac{\partial V}{\partial L}\dot{L}=-[a\left(a+b+e-d_1\right)S^2+\left(a\left(e+c_1\right)+d_1\left(b+e\right)\right)LS+\delta (aS+c_{2}H\\ +\left(a\left(c_1+c_3+d_1-d_3\right)\right)ST+c_2\left(e+c_1\right)HL-c_2\left(a+b-d_1\right)HS\\ -c_2(e-c_3-d_1+d_{3}HT-d_1\left(e-c_1\right)L^2+d_1\left(e-c_1-c_3+d_3\right)LT\end{array}$

$\begin{array}{c}\frac{\partial V}{\partial H}\dot{H}=-[e\left(d_2+c_2-f\right)+d_2+d_2\left(c_1-d_1-d_2\right)HL+e\left(d_1+d_2-c_1-f\right)L^2\\ +\delta eL+e\left(d_2-c_2\right)LS+c_3(f+c_1-d_1-d_2+e\left(c_2-d_2-d_3\right)LT\\ +d_2\left(c_2-d_2-d_3\right)HT+c_3\left(b-d_2\right)ST+c_3\left(d_2-c_2\right)T^2\\ +d_2\left(c_2-f\right)H^2+d_2\left(b-c_2+d_2\right)HS\end{array}$

$\begin{array}{c}\frac{\partial V}{\partial T}\dot{T}=-[f\left(e-c_1-c_3-d_1+d_3\right)HL-f\left(2a+b+c_3+d_3\right)HS+\delta fH\\ +d_3\left(2a+b+c_3-d_3\right)ST+d_3\left(c_3-c_2-d_1-d_3\right)T^2]+d_3(e-c_1-c_3\\ -d_1+d_3)LT+f\left(c_2-c_3-d_2-d_3\right)HT+(\left(2a+b+c_3-d_3\right)S\\ +\delta +\left(c_1+c_3+d_1-d_3\right)L+\left(c_3-c_2-d_1-d_3\right)T)P\left(t,S,L,H,T\right)\end{array}$

$\begin{array}{c}\dot{V}=-[{\tau }_1{S}^2+{\tau }_2{L}^2+{\tau }_3{H}^2+{\tau }_4{T}^2+{\tau }_{5}LS+{\tau }_{6}HL+{\tau }_{7}LT\\ +{\tau }_{8}HS+{\tau }_{9}ST+{\tau }_{10}HT+{\tau }_{11}P\left(t,S,L,H,T\right)\end{array}$

where

${\tau }_1=a\left(a+b\right)+a\left(a+b+e-d_1\right),$

${\tau }_2=\left(a+b-e-d_1\right)c_1+d_1\left(e-c_1\right)+d_1\left(e-c_1\right)+e\left(d_1+d_2-c_1-f\right),$

${\tau }_3=d_2\left(c_2-f\right),$

${\tau }_4=c_3\left(d_2-c_2\right)+d_3\left(c_3-c_2-d_1-d_3\right)$

${\tau }_5={\left(a+b\right)}^2+e\left(a+b\right)+c_1\left(a+b\right)+a\left(e+c_1\right)+d_1\left(b+e\right)+e\left(d_2-c_2\right)$

${\tau }_6=c_1\left(b+c_2-d_3\right)+c_2\left(e+c_1\right)+d_2\left(c_1-d_1-d_2\right)+f\left(e-c_1-c_3-d_1+d_3\right)$

$\begin{array}{c}{\tau }_7=\left(a+b+d_3\right)c_1+d_1\left(e-c_1-c_3+d_3\right)+c_3(f+c_1-d_1-d_2\\ +e\left(c_2-d_2-d_3+d_3\left(e-c_1-c_3-d_1+d_3\right)\right)\end{array}$

$\begin{array}{c}{\tau }_8=c_2\left(a+b\right)+d_2\left(a+b\right)+c_2\left(a+b-d_1\right)+d_2\left(b-c_2+d_2\right)\\ +f\left(2a+b+c_3+d_3\right)\end{array}$

$\begin{array}{c}{\tau }_9={\left(a+b\right)}^2+a\left(a+b\right)+c_3\left(a+b\right)+a\left(c_1+c_3+d_1-d_3\right)+c_3\left(b-d_2\right)\\ +d_3\left(2a+b+c_3-d_3\right)\end{array}$

${\tau }_{10}=c_2\left(e-c_3-d_1+d_3\right)+d_2\left(c_2-d_2-d_3\right)+f\left(c_2-c_3-d_2-d_3\right)$

${\tau }_{11}=\left(2a+b+c_3-d_3\right)S+\delta +\left(c_1+c_3+d_1-d_3\right)L+\left(c_3-c_2-d_1-d_3\right)T$

$\begin{array}{c}\dot{V}=-[{\xi }_1{S}^2+{\xi }_2{L}^2+{\xi }_3{H}^2+{\xi }_4{T}^2+{\tau }_5\left(L^2+S^2\right)+{\tau }_6\left(H^2+L^2\right)\\ +{\tau }_7\left(L^2+T^2\right)+{\tau }_8\left(H^2+S^2\right)+{\tau }_9\left(S^2+T^2\right)+{\tau }_{10}\left(H^2+T^2\right)\\ +{\xi }_5\left(\left|S\right|+\left|L\right|+\left|H\right|+\left|T\right|\right)\mu \left(t\right)\end{array}$

where

$\begin{array}{c}{\xi }_1=\max \{a\left(a+b\right)+a\left(a+b+e-d_1\right),({\left(a+b\right)}^2+e\left(a+b\right)\\ +c_1\left(a+b\right)+a\left(e+c_1\right)+d_1\left(b+e\right)+e\left(d_2-c_2\right),(c_2\left(a+b\right)\\ +d_2\left(a+b\right)+c_2\left(a+b-d_1\right)+d_2\left(b-c_2+d_2\right)\\ +f\left(2a+b+c_3+d_3\right),{\left(a+b\right)}^2+a\left(a+b\right)+c_3\left(a+b\right)\\ +a\left(c_1+c_3+d_1-d_3\right)+c_3\left(b-d_2\right)+d_3\left(2a+b+c_3-d_3\right)\}\end{array}$

$\begin{array}{c}{\xi }_2=\max \{\left(a+b-e-d_1\right)c_1+d_1\left(e-c_1\right)+d_1\left(e-c_1\right)\\ +e\left(d_1+d_2-c_1-f\right),({\left(a+b\right)}^2+e\left(a+b\right)+c_1\left(a+b\right)\\ +a\left(e+c_1\right)+d_1\left(b+e\right)+e\left(d_2-c_2\right),(c_1\left(b+c_2-d_3\right)+c_2\left(e+c_1\right)\\ +d_2\left(c_1-d_1-d_2\right)+f\left(e-c_1-c_3-d_1+d_3\right),(\left(a+b+d_3\right)c_1\\ +d_1\left(e-c_1-c_3+d_3\right)+c_3(f+c_1-d_1-d_2+e(c_2-d_2-d_3\\ +d_3\left(e-c_1-c_3-d_1+d_3\right)\}\end{array}$

$\begin{array}{c}{\xi }_3=\max \{d_2\left(c_2-f\right),(c_1\left(b+c_2-d_3\right)+c_2\left(e+c_1\right)\\ +d_2\left(c_1-d_1-d_2\right)+f\left(e-c_1-c_3-d_1+d_3\right),(c_2\left(a+b\right)\\ +d_2\left(a+b\right)+c_2\left(a+b-d_1\right)+d_2\left(b-c_2+d_2\right)\\ +f\left(2a+b+c_3+d_3\right),c_2\left(e-c_3-d_1+d_3\right)\\ +d_2\left(c_2-d_2-d_3\right)+f\left(c_2-c_3-d_2-d_3\right)\}\end{array}$

$\begin{array}{c}{\xi }_4=\max \{c_3\left(d_2-c_2\right)+d_3\left(c_3-c_2-d_1-d_3\right),\left(a+b+d_3\right)c_1\\ +d_1\left(e-c_1-c_3+d_3\right)+c_3\left(f+c_1-d_1-d_2\right)\\ +e(c_2-d_2-d_3+d_3\left(e-c_1-c_3-d_1+d_3\right),{\left(a+b\right)}^2\\ +a\left(a+b\right)+c_3\left(a+b\right)+a\left(c_1+c_3+d_1-d_3\right)+c_3\left(b-d_2\right)\\ +d_3\left(2a+b+c_3-d_3\right),c_2\left(e-c_3-d_1+d_3\right)\\ +d_2\left(c_2-d_2-d_3\right)+f\left(c_2-c_3-d_2-d_3\right)\}\end{array}$

${\xi }_5=\max \left\{\left(2a+b+c_3-d_3\right),\delta +\left(c_1+c_3+d_1-d_3\right),\left(c_3-c_2-d_1-d_3\right)\right\}$

$\left|P\left(t,S,L,H,T\right)\right|\le \left(\left|S\right|+\left|L\right|+\left|H\right|+\left|T\right|\right)\mu \left(t\right)$

$\dot{V}=-\left[{\xi }_6{S}^2+{\xi }_7{L}^2+{\xi }_8{H}^2+{\xi }_9{T}^2+{\xi }_{10}\left(4+S^2+L^2+H^2+T^2\right)\mu \left(t\right)\right]$

$\dot{V}\le {\xi }_{11}\mu \left(t\right)+{\xi }_{10}\left(S^2+L^2+H^2+T^2\right)\mu \left(t\right)$

integrating the inequality from 0 to t, we have

$V\left(t\right)-V\left(0\right)\le {\xi }_{11}{\displaystyle {\int }_0^t\mu }\left(s\right)\text{d}s+{\xi }_{10}{\displaystyle {\int }_0^{t}V}\left(s\right)\mu \left(s\right)\text{d}s,$

By Gronwall inequality from results of Olutimo and Omoko et al. [15] [20].

$V\left(t\right)\le {\xi }_{12}\exp \left({\xi }_{11}{\displaystyle {\int }_0^{t}V}\left(s\right)\mu \left(s\right)\text{d}s\right)=\text{Γ}$ , Thus the solution of (4) satisfies $\left|S\right|\le \Gamma ,\left|L\right|\le \Gamma ,\left|H\right|\le \Gamma ,\left|T\right|\le \Gamma$ This demonstrates that the total population in each class will stabilize within feasible bounds.

Simulation Results

Wolfram Mathematica provides a powerful symbolic and numerical computing environment for analyzing nonlinear population models giving exact conditions for stability and bounds of solutions. It generate high quality plots for analyzing population trends. Here are some results.

7. Results and Discussions

Construction of the intrinsic Lyapunov functions play a fundamental role in analyzing asymptotic stability and boundedness of population models without solving the equations since they do not have explicit and analytical solutions. This method is applicable to nonlinear systems with several compartments. According to Theorem 1 conditions in the derivation of the Lyapunov function are satisfied, Mathematica’s visualizations in Figure 1 illustrates stability of solutions. Susceptible population gradually decreases as individuals transition into drug use. Light Users, L Peaks and then declines as some progress to heavy use or perhaps, exit the system. Heavy users, H on the other hand, increase initially, stabilize, and decline due to treatment effectiveness T. Figure 2 is the corresponding Wolfram Mathematica code to simulate the stability of the four-compartment drug user model. Figure 3 indicates another population dynamics plot and confirms visual boundedness. The Lyapunov derivative helps mathematically to verify the bounded nature of the solutions. Non-negativity ensures that the populations stay within biologically meaningful bounds. Figure 4 shows the population dynamics of heavy drug users H under different scenarios, Red, Green, Blue, Magenta : The population of heavy users stabilizes at a higher level without any interventions. Increased Treatment and enhancing the treatment rate $c_3$ significantly reduces the steady-state population of heavy users. Reducing the progression rate $d_3$ from lighter to heavy users results in fewer heavy users over time. Figure 5 is a vector stream plot or Phase Plot between light users and heavy users to visualize the population trajectories for stability analysis.

Figure 1. Stability of solutions while increasing $c_3=0.08$ , reducing $d_3=0.06$ , $\delta =100$ .

Figure 2. Wolfram code for stability.

Figure 3. Maintaining the bounds of V(t) at $d_3=0.06,c_3=0.08,\delta =100$ .

Figure 4. Enhancing different drug treatment rate at $c_3=0.1,0.08,0.06,0.04$ .

Figure 5. Stream plot showing interactions between light users, L and heavy users H.

8. Summary & Conclusion

In summary, an intrinsic Lyapunov method defines a Lyapunov stability function that respects the inherent structure of the system, whether geometric, algebraic, or topological. It provides a powerful tool for analyzing complex dynamical systems such as the population models with integration by parts procedure. The dynamics of drug user populations in a heterogeneous population are analyzed using mathematical models that incorporate stability theory. The choice of $V\left(S,H,L,T\right)$ reflects system-specific properties, and its derivative $\dot{V}\left(S,H,L,T\right)$ ensures stability conditions like asymptotic stability to equilibrium. These models often describe how the number of drug users changes over time, influenced by interactions within and between sub populations. Boundedness ensures the model’s solutions remain practically feasible and demonstrate how targeted interventions can stabilize or reduce heavy drug users. This implies that the total population heavy drug users does not grow indefinitely. Stability and Boundedness criteria depend on the balance of parameters such as $c_3,d_3,\delta$ Effective treatment $c_3$ , reduction of relapses $d_3$ , preventive measures and implementing education campaigns targeting youth and vulnerable groups, $\delta$ , all reducing the steady-state $H$ , and consequently, improving public health outcomes. Proper parameter tuning ensures asymptotic stability ($V<0,\dot{V}<0$ ) and minimizes heavy drug users over time. Lyapunov functions evaluate the impact of rehabilitation programs and preventive measures by establishing conditions under which control efforts, such as rehabilitation and awareness programs, effectively eliminate drug use.

Acknowledgements

Sincere thanks to Dr. Olutimo Akinwale Lewis, Lagos State University, Nigeria and special thanks to Mr. Chukwunalu Ajeliba Onwubolu for his immense support always.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Arnold, V.I. (1978) Mathematical Methods of Classical Mechanics, Geometric Mechanics Perspective.
[2]	Hafstein, S. and Giesl, P. (2015) Computational Methods for Lyapunov Functions. Discrete and Continuous Dynamical Systems—Series B, 20, i-ii. https://doi.org/10.3934/dcdsb.2015.20.8i
[3]	Van der Schaft, A.J. (2000) L2-Gain and Passivity Techniques in Nonlinear Control, Application in Control Systems.
[4]	Chin, P.S.M. (1987) Generalized Integral Method to Derive Lyapunov Functions for Nonlinear Systems. International Journal of Control, 46, 933-943. https://doi.org/10.1080/00207178708547404
[5]	Lin, S.-Z., Liu, Z.-R. and Yu, Y.-H. (1998) Stability for Certain Fourth Order Nonlinear Differential Equations. Demonstratio Mathematica, 31, 87-96.
[6]	Chin, P.S.M. (1988) Extension to the Intrinsic Method and Its Application to Stability Problems. International Journal of Control, 48, 2139-2146. https://doi.org/10.1080/00207178808906312
[7]	Chin, P.S.M. (1988) Stability of Non-Linear Systems via the Intrinsic Method. International Journal of Control, 48, 1561-1567. https://doi.org/10.1080/00207178808906269
[8]	Tunç, C. (2007) Some Remarks on the Stability and Boundedness of Solutions of Certain Differential Equations of Fourth-Order. Computational & Applied Mathematics, 26, 1-17. https://doi.org/10.1590/s0101-82052007000100001
[9]	Lucas Matonya, F. and Kuznetsov, D. (2019) Mathematical Modelling on the Effects of Drug Abuse to the Societies. Annals of Pure and Applied Mathematics, 19, 21-35. https://doi.org/10.22457/apam.591v19n1a4
[10]	Peranginangin, A.P. (2023) Stability Analysis of Nonlinear Dynamic Systems Using the Lyapunov Method Approach. Mandalika Mathematics and Educations Journal, 5, 253-265. https://doi.org/10.29303/jm.v5i2.6383
[11]	Matonya, F.L. and Kuznetsov, D. (2021) Mathematical Modelling of Drug Abuse and It’s Effect in the Society. Annals of Pure and Applied Mathematics, 24, 159-176. https://doi.org/10.22457/apam.v24n2a06855
[12]	Resmawan, R., Supu, N. and Achmad, N. (2021) Dynamic Model Analysis of the Spread of Drug Addicts with Educational Effects. In: Proceedings of the 1st International Conference on Mathematics and Mathematics Education, Atlantis Press, 55-64. https://doi.org/10.2991/assehr.k.210508.042
[13]	Oni, O. and Awofala, T.B. (2022) Utilization of Wearable Technology Data in Chronic Disease Management. World Journal of Advanced Research and Reviews, 16, 1189-1195. https://doi.org/10.30574/wjarr.2022.16.3.1208
[14]	Olutimo, A.L., Oni, O.J., Williams, F.A., Akewushola, J.R. and Abass, F.A. (2024) Transmission and Control Dynamics of Rotavirus Diarrhea Model with Double Dose Vaccination. Journal of Mechanics of Continua and Mathematical Sciences, 19, 67-86. https://doi.org/10.26782/jmcms.2024.10.00005
[15]	Omoko, I.D., Abdurasid, A.A., Akewushola, J.R. and Salami, O.A. (2024) Qualitative Analysis of Solutions of a Nonlinear Biomechanical Model. Journal of Applied Mathematics and Physics, 12, 3982-3993. https://doi.org/10.4236/jamp.2024.1211242
[16]	Ogundare, B.S. (2009) Qualitative and Quantitative Properties of Solutions of Ordinary Differential Equations. PhD Thesis, Department of Pure and Applied Mathematics, University of Fort Hare Alice.
[17]	Lyapunov, A.M. (1992) The General Problem of the Stability of Motion. Translated and Edited by A.T. Fuller, Taylor and Franscis.
[18]	LaSalle, J.P. (1968) The Stability of Dynamical Systems.
[19]	Olutimo, L.A., Oni, O.J. and Omoko, I.D. (2020) On the Boundedness Analysis of the State Variables for a Loaded DC Servo Motor System. 2020 International Conference in Mathematics, Computer Engineering and Computer Science (ICMCECS), Ayobo, 18-21 March 2020, 1-4. https://doi.org/10.1109/icmcecs47690.2020.240899
[20]	Olutimo, A.L. and Omoko, I.D. (2020) Stability and Boundedness Analysis of a System of RLC Circuit with Response. International Journal of Applied Mathematics, 33, 39-48. https://doi.org/10.12732/ijam.v33i1.4