Mathematical modeling of carbon dioxide emissions with GDP linkage: sensitivity analysis and optimal control strategy

TL;DR

This paper uses a mathematical model to study how GDP, forests, and population affect CO₂ levels, finding that forests have a bigger impact than GDP on reducing emissions.

Contribution

The novelty lies in using a mathematical model with sensitivity analysis to compare the influence of GDP and forests on CO₂ emissions.

Findings

GDP-driven CO₂ emissions and atmospheric CO₂ concentration have limited sensitivity in the model.

Forests show higher sensitivity and greater influence on CO₂ levels compared to GDP.

Optimal control strategies suggest prioritizing forestation for effective CO₂ mitigation.

Abstract

Climate change and global warming are among the most significant issues that humanity is currently facing, and also among the issues that pose the greatest threats to all mankind. These issues are primarily driven by abnormal increases in greenhouse gas concentrations. Mathematical modeling serves as a powerful approach to analyze the dynamic patterns of atmospheric carbon dioxide. In this paper, we established a mathematical model with four state variables to investigate the dynamic behavior of the interaction between atmospheric carbon dioxide, GDP, forest area and human population. Relevant theories were employed to analyze the system’s boundedness and the stability of equilibrium points. The parameter values were estimated with the help of the actual data in China and numerical fitting was carried out to verify the results of the theoretical analysis. The Partial Rank Correlation Coefficient (PRCC) determines the sensitivity ofan input parameter to the output by measuring the correlation between a single input parameter and the model output. The sensitivity analysis of the compartments with respect to the model parameters was analyzed by using the PRCCand the Latin Hypercube Sampling test.The results indicate that the sensitivity of GDP-driven CO₂ emissions and GDP-governed atmospheric CO₂ concentration to the system is not significant. This implies that within the GDP-driven mitigation framework, the regulatory effect of GDP on atmospheric CO₂ concentration is relatively limited, and its significance is less pronounced than that of forests. Therefore, future relevant strategies should prioritize parameters with higher sensitivity (e.g., forestation). Apply the optimal control theory to regulate the atmospheric carbon dioxide level and provide the corresponding numerical fitting. Finally, corresponding discussions and suggestions were put forward with the help of the results of the theoretical analysis and numerical fitting.

Figures

Trends of CO_2_ [[5](#CR5)]

Flow chart of the model system (5)

Global Stability of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_4}$$\end{document}$ in *C-G-F* space

Global Stability of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_4}$$\end{document}$ in *N-F-C* space

Time series graph of *C* and *F* for different value of parameter $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document}$

Time series graph of *C* and *N* for different value of parameter $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document}$

The PRCC value of each parameter with respect to compartment C(t), G(t), F(t) and N(t) at t = 4000 respectively

The variation of PRCC value of the parameter with respect to (a) carbon dioxide, (b) GDP, (c) forest area and (d) human population

Graph trajectories with and without optimal control for: (a) carbon dioxide, (b) GDP, (c) forest area, (d) human population

Full text

Introduction

Since the commencement of the Industrial Revolution, the accelerated pace of industrialization in human society has notably aggravated the emission of greenhouse gases, with carbon dioxide (CO₂) emerging as the principal contributor to this phenomenon. This acceleration, driven by exponential growth in fossil fuel combustion, industrial processes, and deforestation, has disrupted the natural carbon cycle, thereby exacerbating the greenhouse effect and altering global climatic patterns [1]. Throughout this epoch, anthropogenic activities, principally fossil fuel combustion, deforestation, and land-use transformation have propelled atmospheric CO₂ concentrations from approximately 280 parts per million (ppm) to over 400 ppm. This unprecedented rise, equivalent to a 43% increase since the pre-industrial era, reflects a systematic disruption of the global carbon cycle, with profound implications for climate stability and ecological balance [2]. This increase in concentrations is the main driver of global warming, which exacerbates the rise in global temperatures and triggers widespread climate change phenomena [3]. By February 2025, the atmospheric carbon dioxide concentration had surged to an astonishing 426.13 parts per million (ppm) [4]. Fig. 1 depicts the temporal dynamics of annual global greenhouse gas emissions from 1980 to 2023, revealing a consistent upward trajectory. While inter-annual fluctuations are evident, the secular growth trend remains statistically robust, underscoring the escalating severity of greenhouse gas emissions and the urgent imperatives for climate change mitigation. This trend not only reflects the cumulative impact of anthropogenic activities but also highlights the critical need for concerted, systemic action to address this planetary challenge [5].Fig. 1. Trends of CO_2_ [5]

Climate change is manifested in rising global average temperatures, frequent extreme weather events, accelerating glacier melting, continuous sea-level rise, and significant ecosystem disruptions [6]. These changes pose a serious threat to global agricultural production, water allocation, biodiversity conservation and human health [7]. For example, sea-level rise endangers coastal communities and infrastructure, extreme weather intensifies natural disaster risks, and ecosystem disruptions alter species distributions and survival rates [8].

In response to global warming and climate change, the international community has taken various actions. The Paris Agreement stands as a key milestone, aiming to keep global temperature rise well below 2 °C and pursue efforts for 1.5 °C [9]. Countries are adopting multiple strategies to cut greenhouse gas emissions: developing renewable energy, boosting energy efficiency, implementing emissions trading schemes, and restoring forests [10]. The alterations in extreme weather conditions, coupled with the rise in global surface temperature and climate change, have given rise to the spread of a host of diseases [1]. Climate change is evidenced by multiple phenomena: an increase in global mean temperature, a surge in extreme weather occurrences, glacial and ice sheet melt, sea-level ascent, and significant ecosystem alterations [11]. Mathematical modeling can effectively visualize the dynamic behavior of atmospheric carbon dioxide, enabling better adoption of corresponding measures to alleviate the level of atmospheric carbon dioxide [12].

Numerous mathematical frameworks have been put forward to investigate the impacts of factors such as carbon capture technology, population pressure, reforestation, vehicle CO_2_ emissions, technology choice and urbanization on the dynamics of atmospheric carbon dioxide [13–23]. In particular, Devi and Gupta [24] proposed a nonlinear mathematical model to simulate changes in the ability of plants to absorb atmospheric carbon dioxide. The paper indicates that afforestation represents a scientific approach to lower atmospheric carbon dioxide concentrations. In [25], a three-dimensional mathematical model has been constructed to analyze the impact of budget al.location on the reduction of atmospheric carbon dioxide concentrations. As the level of carbon dioxide in the atmosphere increases, the growth rate of budget al.location may lead to a stability switch through the Hopf-bifurcation. Mishra et al. [26] discussed the use of green belt planting and seaweed farming to reduce atmospheric carbon dioxide (CO_2_), and model analysis showed that the use of plants for photosynthesis by planting leafy trees in the green belt around the emission source, and through seaweed farming, could effectively reduce atmospheric CO_2_ levels. Tandon [27] employed a mathematical model to investigate the impacts of mining activities on the dynamic natural interactions between plants and carbon dioxide. The study revealed that mining activities notably elevated atmospheric carbon dioxide concentrations and caused damage to plants, thereby hindering the system’s ability to attain a stable state. Misra and Verma [28] studied the effects of population and forest biomass on atmospheric carbon dioxide. The results show that when human deforestation exceeds a certain threshold, the system will occur Hopf-bifurcation and become unstable. Most of the above-mentioned literatures have considered the kinetic relationship between carbon dioxide and forest. Caetano et al. [29] focused their analysis on how GDP affects atmospheric carbon dioxide, neglecting to consider the interrelationship between carbon dioxide and forest ecosystems. Our modeling framework accounts for the interactive relationship between CO₂ emissions and GDP, constructing a four-dimensional mathematical model to capture the dynamic behavior of atmospheric carbon dioxide.

China is the country with the largest carbon emissions in the world, accounting for nearly one-third of the global total carbon emissions. And this is inextricably linked to China’s rapidly growing economy. Therefore, researching the dynamic relationship of the interaction between China’s GDP index and the carbon cycle is of far-reaching significance in the efforts to mitigate global warming and climate change.

Mathematical model

In this work, we proposed a mathematical model to research the dynamics of CO_2_ emissions, GDP, forest area and human population. The variables are defined as follows:

1. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C\left( t \right)$$\end{document} : the concentration of atmospheric CO_2_ (in ppm).
2. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G\left( t \right)$$\end{document} : the gross domestic product (in billion USD$).
3. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\left( t \right)$$\end{document} : the forest area (in million hectares.).
4. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\left( t \right)$$\end{document} : the human population (in million).

The atmospheric carbon dioxide emissions stem from two categories: The emission of natural factors (such as volcanic eruptions, respiration processes of plants and animals, etc.), this constant growth term we denote as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} . The emission of carbon dioxide caused by human activities which is proportional to the population [30], we denote this increase as parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} . Rapid economic growth is often accompanied by a lot of industrialization [31], i.e. more carbon dioxide is emitted into the atmosphere, which we record as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} . Forest area absorbs atmospheric carbon dioxide through the process of photosynthesis and leads a decrease in atmospheric carbon dioxide [32], which we denote as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} . The lifetime of atmospheric carbon dioxide is usually 30 to 95 years [33], we denote the natural loss coefficient of atmospheric carbon dioxide as p. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} represent the growth rate of GDP. Economic growth, in turn, can reduce carbon dioxide in the atmosphere through activities such as clean technologies, we use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document} to represent it. Based on these assumptions, the dynamics of atmospheric carbon dioxide are governed by the following equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dC}}{{dt}}=\alpha +\phi N{\text{+}}\left( {\beta - \varepsilon } \right)G - \eta CF - p C$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dG}}{{dt}}=\mu - \varepsilon G$$\end{document}

In the modeling process, we assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} represent the natural growth rate of forest areas and K represent the carrying capacity of forest area. Human population growth often leads to an increase in demand for deforestation, which may include agricultural land, urban sprawl and infrastructure development. To meet these needs, forests may be cut down and leading to a decrease in forest area [34], this anthropogenic deforestation coefficient is denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} . The absorption of the right amount of carbon dioxide will promote the growth of forest area more densely [35]. The growth rate of forest area caused by the absorption of carbon dioxide we denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} . Based on these assumptions, the dynamics of forest area are governed by the following equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dF}}{{dt}}=\omega F\left( {1 - \frac{F}{K}} \right) - \theta NF+\eta \sigma CF$$\end{document}

We assume that the population follows logistic growth, with s and M denoting the natural growth rate of the population and the carrying capacity of the population respectively. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} represents the contribution of forest area to the population (e.g. provision of food and resources, conditions, climate, etc.) [18]. The absorption of carbon dioxide by the human body can have dire consequences, it can directly lead to death or exacerbate the spread of specific diseases, we use parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} to express the rate of loss of the population caused by carbon dioxide [23]. Based on these assumptions, the dynamics of popultion governed by the following equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dN}}{{dt}}=sN\left( {1 - \frac{N}{M}} \right)+\theta \nu NF - \pi CN$$\end{document}

In summary, our mathematical model is as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \frac{{dC}}{{dt}}=\alpha +\phi N{\text{+}}\left( {\beta - \varepsilon } \right)G - \eta CF - pC \hfill \\ & \frac{{dG}}{{dt}}=\mu - \varepsilon G \hfill \\ & \frac{{dF}}{{dt}}=\omega F\left( {1 - \frac{F}{K}} \right) - \theta NF+\eta \sigma CF \hfill \\ & \frac{{dN}}{{dt}}=sN\left( {1 - \frac{N}{M}} \right)+\theta \nu NF - \pi CN \hfill \\ \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C\left( 0 \right)={C_0} \geq 0$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G\left( 0 \right)={G_0} \geq 0$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\left( 0 \right)={F_0} \geq 0$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\left( 0 \right)={N_0} \geq 0$$\end{document} . Fig. 2 shows the flow chart of the model system (5).Fig. 2. Flow chart of the model system (5)

Model analysis

Boundedness

The boundedness of the system is given by the following lemma:

Theorem 3.1.1

If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha +\phi {N_{max}}$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+\left( {\beta - \varepsilon } \right)$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G_{max}}>0$$\end{document} , then solutions of system (5) are bounded in region  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega =\left\{ \left( {C,G,F,N} \right) \in R_{+}^{4}: \right.$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left. 0 \leq C \leq {C_{max}};\right.$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left. 0 \leq G \leq {G_{max}};\right.$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left. 0 \leq F \leq {F_{max}}; \right.$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left. 0 \leq N \leq {N_{max}} \right\}$$\end{document} , here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_{\hbox{max} }}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G_{\hbox{max} }}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${F_{\hbox{max} }}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${N_{\hbox{max} }}$$\end{document} are given as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_{max}}=\frac{{\alpha +\phi {N_{max}}+\left( {\beta - \varepsilon } \right){G_{max}}}}{p}$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G_{max}}=\frac{\mu }{\varepsilon }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${F_{max}}=\frac{{K\left( {\omega +\eta \sigma {C_{max}}} \right)}}{\omega }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${N_{max}}=M+\frac{{\theta \nu M}}{s}{F_{max}}$$\end{document} and this attracts all solutions initiating from within the positive orthant’s interior.

Proof According to the comparison theorem we get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 \leq G\left( t \right) \leq \frac{\mu }{\varepsilon }={G_m}_{{{\text{ax}}}}(say)$$\end{document}

From the first equation we get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dC}}{{dt}} \leq \alpha +\phi N+\left( {\beta - u} \right){G_{{max} }} - pC$$\end{document}

which gives that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathop {\lim \sup C(t)}\limits_{{t \to \infty }} & =\frac{{\alpha +\phi {N_{max}}+\left( {\beta - u} \right){G_{max}}}}{p} \\ & ={C_m}_{{{\text{ax}}}}(say) \end{aligned}$$\end{document}

From the third equation we get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dF}}{{dt}} \leq \omega F\left( {1 - \frac{F}{K}} \right)+\eta \sigma CF+\varepsilon G$$\end{document}

which implies that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {\lim \sup F(t)}\limits_{{t \to \infty }} =\frac{{K\left( {\omega +\eta \sigma {C_{max}}} \right)}}{\omega }={F_m}_{{{\text{ax}}}}(say)$$\end{document}

From the fourth equation we get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dN}}{{dt}} \leq \left( {s+\theta \nu F} \right)N - \frac{s}{M}{N^2}$$\end{document}

which implies that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathop {\lim \sup N(t)}\limits_{{t \to \infty }} & =\frac{{M(s+\theta \nu {F_{max}})}}{s}\\ & =M+\frac{{\theta \nu }}{s}{F_{max}}={N_m}_{{{\text{ax}}}}(say) \end{aligned}$$\end{document}

This completes the proof of the boundedness of the system (5).

Equilibrium points

To solve the equilibrium point of the model, we need to make the right of the model equal to zero to solve all possible equilibrium points. The model system (5) has two non-negative equilibria which are listed as follows:

• (i) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{1}}}=\left( {{C_1},{G_1},0,0} \right)$$\end{document}
• (ii) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{2}}}=\left( {{C_{\text{2}}},{G_{\text{2}}},{F_{\text{2}}},0} \right)$$\end{document}
• (iii) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{3}}}=\left( {{C_{\text{3}}},{G_{\text{3}}},0,{N_{\text{3}}}} \right)$$\end{document}
• (iv) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{4}}}=\left( {{C_{\text{4}}},{G_{\text{4}}},{F_{\text{4}}},{N_{\text{4}}}} \right)$$\end{document}

Theorem 3.2.1

System (5) possesses an equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_1}$$\end{document} provided that the following inequality is satisfied:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha {\text{+}}\left( {1 - \frac{\beta }{\varepsilon }} \right)>0$$\end{document}

Proof

• (i)From equation two we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dG}}{{dt}}=\mu - \varepsilon G=0$$\end{document}

therefor we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G_{\text{1}}}{\text{=}}{G_*}=\frac{\mu }{\varepsilon }$$\end{document}

When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=F=0$$\end{document} , from equation one we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C=\frac{1}{p}\left[ {\alpha +\left( {\frac{\beta }{\varepsilon } - \mu } \right)} \right]$$\end{document}

Thus we get equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_1}=\left( \frac{1}{p}\left[ \vphantom{\frac{\beta }{\varepsilon }} \alpha \right.\right.$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left.\left. +\left( \frac{\beta }{\varepsilon } - \mu \right) \right], \frac{\mu }{\varepsilon },0,0 \right)$$\end{document}

Theorem 3.2.2

System (5) possesses an equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_2}$$\end{document} provided that the following inequality is satisfied:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \varepsilon +\mu \left( {\beta - \varepsilon } \right)>0$$\end{document}

• (ii)When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N={\text{0}}$$\end{document} , from equation one and three we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega -\frac{\omega }{K}F+\eta \sigma C=0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha +\left( {\beta - \varepsilon } \right)\frac{\mu }{\varepsilon } - \eta CF - pC=0$$\end{document}

Using the value of F from (11) in (12), we get the following quadratic polynomial in C

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{{\eta ^2}\sigma K}}{\omega }{C^2}+\left( {\eta K+p} \right)C - \left[ {\alpha - \left( {1 - \frac{\beta }{\varepsilon }} \right)} \right]=0 \quad$$\end{document}

Hence, applying Descartes’ rule of signs confirms the existence of a unique positive root

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_2}=\frac{{\omega \left[ { - \left( {p+\eta K} \right)+\sqrt {\left( {{{\left( {p+\eta K} \right)}^2}+\frac{{4\sigma {\eta ^2}K}}{\omega }\left[ {\alpha +\left( {1 - \frac{\beta }{\varepsilon }} \right)} \right]} \right)} } \right]}}{{2\sigma {\eta ^2}K}} \quad$$\end{document}

if condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \varepsilon +\mu \left( {\beta - \varepsilon } \right)>0$$\end{document} satisfied. Substitute Eq. (14) into Eq. (11), and a unique positive value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${F_2}$$\end{document} can be obtained.

Theorem 3.2.3

System (5) possesses an equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_3}$$\end{document} provided that the following inequality is satisfied:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon \left( {\alpha +\phi M} \right)+\mu \left( {\beta - \varepsilon } \right)>0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon \left( {sp+\pi \phi M} \right)>\varepsilon \left( {\alpha +\phi M} \right)+\mu \left( {\beta - \varepsilon } \right)$$\end{document}

• (iii)when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0$$\end{document} , from equation one and four of systsem (5) we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha +\phi N{\text{+}}\left( {\beta - \varepsilon } \right)G - pC=0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s - \frac{s}{M}N - \pi C=0$$\end{document}

Using the value of N from (18) in (17) we get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C=\frac{{s\left[ {\varepsilon \left( {\alpha +\phi M} \right)+\mu \left( {\beta - \varepsilon } \right)} \right]}}{{\varepsilon \left( {sp+\pi \phi M} \right)}}$$\end{document}

Using the value of C in (18) we get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=M\left\{ {1 - \frac{{\pi \left[ {\varepsilon \left( {\alpha +\phi M} \right)+\mu \left( {\beta - \varepsilon } \right)} \right]}}{{\varepsilon \left( {sp+\pi \phi M} \right)}}} \right\}$$\end{document}

Theorem 3.2.4

System (5) possesses an equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_4}$$\end{document} provided that the following inequality is satisfied:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \left[ {\varepsilon p\left( {sp+\pi \phi M} \right)} \right]\left\{ {\varepsilon \omega p+\eta \sigma \left[ {\varepsilon \alpha +\mu \left( {\beta - \varepsilon } \right)} \right]} \right\} \\ & \qquad>\varepsilon M\left( {p\theta - \eta \sigma \phi } \right)\left[ {\varepsilon \left( {\alpha +s\pi p} \right)+\mu \left( {\beta - \varepsilon } \right)} \right] \end{aligned} \quad$$\end{document}

(iv) From system (5) we get following Eq.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha +\phi N+\left( {\beta - \varepsilon } \right)\frac{\mu }{\varepsilon } - \eta CF - pC=0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega - \frac{\omega }{K}F - \theta N+\eta \sigma C=0$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s - \frac{s}{M}N+\theta \nu F - \pi C=0$$\end{document}

From (22) we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C=\frac{{\alpha +\phi N+\left( {\beta - \varepsilon } \right)\frac{\mu }{\varepsilon }}}{{p+\eta F}}$$\end{document}

Using (25) in (23) and (24) we get following Eq.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a:\omega - \frac{\omega }{K}F - \theta N+\eta \sigma \frac{{\alpha +\phi N+\left( {\beta - \varepsilon } \right)\frac{\mu }{\varepsilon }}}{{p+\eta F}}=0 \quad$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b:s - \frac{s}{M}N+\theta \nu F - \frac{{\alpha +\phi N+\left( {\beta - \varepsilon } \right)\frac{\mu }{\varepsilon }}}{{p+\eta F}}=0 \quad$$\end{document}

To prove the existence of equilibrium points, we now analyze curves a and b separately.

For Eq. (26):

• (i)when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F={F_{\text{2}}}>0$$\end{document} if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \varepsilon +\mu \left( {\beta - \varepsilon } \right)>0.$$\end{document} .
• (ii)when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0$$\end{document} , we have \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N={N_a}=$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{\varepsilon \omega p+\eta \sigma \left[ {\varepsilon \alpha +\mu \left( {\beta - \varepsilon } \right)} \right]}}{{\varepsilon \left( {p\theta - \eta \sigma \phi } \right)}}>0$$\end{document} if

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p\theta - \eta \sigma \phi>0$$\end{document}

• (iii)By differentiating Eq. (26) with respect to F, we obtain:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{{dN}}{{dF}} & =\frac{{\varepsilon \omega {{\left( {p+\eta F} \right)}^2}+\sigma {\eta ^2}\left[ {\varepsilon \left( {\alpha +\phi N} \right)+\mu \left( {\beta - \varepsilon } \right)} \right]K}}{{\varepsilon K\left( {p+\eta F} \right)\left[ {\eta \sigma \phi - \theta \left( {p+\eta F} \right)} \right]}} \quad \\ & <0 \end{aligned}$$\end{document}

For Eq. (27):

• (i)when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F=0$$\end{document} , we have \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N={N_b}=$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{M\left[ {\varepsilon \left( {\alpha +s\pi p} \right)+\mu \left( {\beta - \varepsilon } \right)} \right]}}{{\varepsilon p\left( {sp+\pi \phi M} \right)}}>0$$\end{document} if

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon \left( {\alpha +s\pi p} \right)+\mu \left( {\beta - \varepsilon } \right)>0$$\end{document}

• (ii)when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=0$$\end{document} , we get following equation in F

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \nu \varepsilon \eta \theta {F^2}+s\varepsilon \left( {\eta +\nu p} \right)F \\ & \quad +\left[ {\varepsilon \left( {sp - \alpha } \right) - \left( {\beta - \varepsilon } \right)} \right]=0 \end{aligned}$$\end{document}

applying Descartes’ rule of signs confirms the existence of a unique negative root \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F={F_b}<0$$\end{document} if

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon \left( {sp - \alpha } \right) - \left( {\beta - \varepsilon } \right)<0$$\end{document}

• (iii)Calculate the derivative of (27) with respect to F we get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{{dN}}{{dF}} & =\frac{{\nu \varepsilon \theta {{\left( {p+\eta F} \right)}^2}+\eta \left[ {\varepsilon \left( {\alpha +\phi N} \right)+\mu \left( {\beta - \varepsilon } \right)} \right]}}{{\varepsilon \left( {p+\eta F} \right)\left[ {\phi M+s\left( {p+\eta F} \right)} \right]}}M \quad \\ &>0 \end{aligned}$$\end{document}

This indicates that there is a unique intersection between the two curves at equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{4}}}$$\end{document} , at this time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${N_a}>{N_b}$$\end{document} must be satisfied, i.e.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \left[ {\varepsilon p\left( {sp+\pi \phi M} \right)} \right]\left\{ {\varepsilon \omega p+\eta \sigma \left[ {\varepsilon \alpha +\mu \left( {\beta - \varepsilon } \right)} \right]} \right\} \\ & \qquad>\varepsilon M\left( {p\theta - \eta \sigma \phi } \right)\left[ {\varepsilon \left( {\alpha +s\pi p} \right)+\mu \left( {\beta - \varepsilon } \right)} \right] \end{aligned}$$\end{document}

Stability analysis

We discuss the stability of equilibria \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_1}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_2}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_3}$$\end{document} by finding the sign of the eigenvalues of Jacobian matrix corresponding to each equilibrium. The Jacobian matrix for model system (5) is given as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=\left[ {\begin{array}{*{20}{c}} { - \eta F - p}&{\beta - \varepsilon }&{ - \eta C}&\phi \\ 0&{ - \varepsilon }&0&0 \\ {\eta \sigma F}&0&{\omega - \frac{{2\omega }}{K}F - \theta N+\eta \sigma C}&{ - \theta F} \\ { - \pi N}&0&{\nu \theta N}&{s - \frac{{2s}}{M}N+\nu \theta F - \pi C} \end{array}} \right]$$\end{document}

Define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J_i}(i=1,2,3,4)$$\end{document} as the equilibrium Jacobian.

Theorem 4.1

• (i) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_1}$$\end{document} is inherently unstable under all conditions.
• (ii) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_2}$$\end{document} is always stable in G direction and locally stable (unstable) manifold in N direction provided \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s+\nu \theta {F_2} - \pi {C_2}$$\end{document} is negative (positive). Asymptotically stable in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C - F$$\end{document} directions when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega - \frac{{2\omega }}{K}{F_2}+\eta \sigma {C_2}<\hbox{min} \left\{ {\eta {F_2}+p,\frac{{\sigma {\eta ^2}{F_2}{C_2}}}{{\eta {F_2}+p}}} \right\}$$\end{document}
• (iii) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_3}$$\end{document} is always stable in G direction, whereas \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_3}$$\end{document} is locally stable (unstable) manifold in F direction provided \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega - \frac{{2\omega }}{K}{F_2}+\eta \sigma {C_2}$$\end{document} is negative (positive). asymptotically stable in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C - F$$\end{document} directions when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s - \frac{{2s}}{M} - \pi {C_3}<\hbox{min} \left\{ {p,\frac{{\pi \phi {N_3}}}{p}} \right\}$$\end{document}

Proof

• (i)The eigenvalues of the Jacobian matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J_1}$$\end{document} are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- p$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- \varepsilon$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega +\eta \sigma \left[ {\frac{1}{p}+\left( {\frac{\beta }{\varepsilon } - \mu } \right)} \right]$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- \pi \left[ {\frac{1}{p}+\left( {\frac{\beta }{\varepsilon } - \mu } \right)} \right]$$\end{document}  seperately. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega +\eta \sigma \left[ {\frac{1}{p}+\left( {\frac{\beta }{\varepsilon } - \mu } \right)} \right]>0$$\end{document} whenever \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_1}$$\end{document} exist.
• (ii)The eigenvalues of the Jacobian matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J_2}$$\end{document} in G and N directions are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- \varepsilon$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s+\nu \theta {F_2} - \pi {C_2}$$\end{document} , therefore \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_2}$$\end{document} is always stable in G direction, whereas locally stable (unstable) manifold in N direction provideds \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s+\nu \theta {F_2} - \pi {C_2}$$\end{document} is negative (positive). The other two eigenvalues are solutions to unary quadratic equation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${y^2} - \left( \omega +\eta \sigma {C_2} - \eta {F_2} - p - \frac{{2\omega }}{K}{F_2} \right)+\left[ \sigma {\eta ^2}{F_2}{C_2} - \vphantom{\frac{1}{2}} \right.$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left. \left( \eta {F_2}+p \right)\left( \omega - \frac{{2\omega }}{K}{F_2}+\eta \sigma {C_2} \right) \right]=0$$\end{document} . When the trace is less than zero and the value of the determinant is greater than zero, there are negative eigenroots or negative real part, i.e.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega - \frac{{2\omega }}{K}{F_2}+\eta \sigma {C_2}<\hbox{min} \left\{ {\eta {F_2}+p,\frac{{\sigma {\eta ^2}{F_2}{C_2}}}{{\eta {F_2}+p}}} \right\}$$\end{document}

Thus the eigenroot is asymptotically stable when (35) satisfied.

• (iii)The eigenvalues of the Jacobian matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J_3}$$\end{document} in G and F directions are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- \varepsilon$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega - \theta {N_2}+\eta \sigma {C_2}$$\end{document} , therefore \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_3}$$\end{document} is always stable in G direction whereas \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_3}$$\end{document} is locally stable (unstable) manifold in F direction provideds \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega - \frac{{2\omega }}{K}{F_2}+\eta \sigma {C_2}$$\end{document} is negative (positive). The other two eigenvalues are solutions to unary quadratic equation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${y^2} -$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( s - p - \frac{{2s}}{M} - \pi {C_3} \right)+\left[ \pi \phi {N_3} - p\left( s - \frac{{2s}}{M} - \pi {C_3} \right) \right]$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=0$$\end{document} . When the trace is less than zero and the value of the determinant is greater than zero, there are negative eigenroots or negative real part, i.e.
• (iv)

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s - \frac{{2s}}{M} - \pi {C_3}<\hbox{min} \left\{ {p,\frac{{\pi \phi {N_3}}}{p}} \right\}$$\end{document}

Thus the eigenroot is asymptotically stable when (36) satisfied.

Theorem 4.2

Local asymptotic stability of the interior equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{4}}}$$\end{document} is guaranteed if the subsequent condition holds:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A_1}{A_2} - {A_3}>0$$\end{document}

• (iv)**Proof **Evaluating the Jacobian matrix at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{4}}}$$\end{document} we derive:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J_{\text{4}}}=\left[ {\begin{array}{*{20}{c}} { - \eta {F_{\text{4}}} - p}&{\beta - \varepsilon }&{ - \eta {C_{\text{4}}}}&\phi \\ 0&{ - \varepsilon }&0&0 \\ {\eta \sigma {F_{\text{4}}}}&0&{\omega - \frac{{2\omega }}{K}{F_{\text{4}}} - \theta {N_{\text{4}}}+\eta \sigma {C_{\text{4}}}}&{ - \theta {F_{\text{4}}}} \\ { - \pi {N_{\text{4}}}}&0&{\nu \theta {N_{\text{4}}}}&{s - \frac{{2s}}{M}{N_{\text{4}}}+\nu \theta {F_{\text{4}}} - \pi {C_{\text{4}}}} \end{array}} \right]$$\end{document}

The characteristic equation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J_{\text{4}}}$$\end{document} is.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( {\psi +\varepsilon } \right)\left( {{\psi ^3}+{A_1}{\psi ^2}+{A_2}\psi +{A_3}} \right)=0$$\end{document}

where,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A_1}=\eta {F_{\text{4}}}+p+\frac{{2\omega }}{K}{F_{\text{4}}}+2\nu \theta {F_{\text{4}}}+\frac{s}{M}{N_{\text{4}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A_2}=\frac{\omega }{K}{F_{\text{4}}} \cdot \frac{s}{M}{N_{\text{4}}}+\left( {p+\eta {F_{\text{4}}}} \right)\left( {\frac{\omega }{K}{F_{\text{4}}}+\frac{s}{M}{N_{\text{4}}}} \right)$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A_3}=\frac{s}{M}\eta \theta \left( {\pi {C_{\text{4}}}+\nu \sigma \phi } \right){N_{\text{4}}}^{2}{F_{\text{4}}}+\left[ \frac{\omega }{K}\pi \phi +\nu {\theta ^2}\left( \eta {F_{\text{4}}} \right. \right.$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left. \left. +p \right) \right]{N_{\text{4}}}{F_{\text{4}}}+\frac{\omega }{K} \cdot \frac{s}{M}{\eta ^2}\sigma {N_{\text{4}}}{F_{\text{4}}}^{2}$$\end{document}

Here, it is straightforward to see that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A_1}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A_2}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A_3}$$\end{document} are positive. The Routh–Hurwitz criterion therefore simplifies to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A_1}{A_2} - {A_3}>0$$\end{document} .

Theorem 4.3

Global stability of the interior equilibrium \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{4}}}$$\end{document} inside the region of attraction is guaranteed if the subsequent conditions hold:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox{max} \left\{ {\frac{{\eta K}}{{{m_1}\omega }}{{\left( {{m_1}\sigma - {C_{\hbox{max} }}} \right)}^2},\frac{{{{\left( {\beta - \varepsilon } \right)}^2}}}{\varepsilon }} \right\}<2\left( {p+\eta {F_{\text{4}}}} \right)$$\end{document}

Proof To establish the global stability of interior equilibrium \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{4}}}$$\end{document} , we employ Lyapunov’s method by selecting a positive definite function as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}V & =\frac{1}{2}{\left( {C - {C_{\text{4}}}} \right)^2}+\frac{1}{2}{\left( {G - {G_{\text{4}}}} \right)^2}\\ & \quad +{m_1}\left( {F - {F_{\text{4}}} - {F_{\text{4}}}ln\frac{F}{{{F_{\text{4}}}}}} \right)\\ & \quad +{m_2}\left( {N - {N_{\text{4}}} - {N_{\text{4}}}ln\frac{N}{{{N_{\text{4}}}}}} \right) \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${m_1}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${m_2}$$\end{document} represent positive constants subject to appropriate calibration.

Differentiating V with respect to t along the solution path of model system (5) yields:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{{dV}}{{dt}} & = - \left( {p+\eta {F_{\text{4}}}} \right){\left( {C - {C_{\text{4}}}} \right)^2} - \varepsilon {\left( {G - {G_{\text{4}}}} \right)^2} \\ & \quad - \frac{\omega }{K}{\left( {F - {F_{\text{4}}}} \right)^2} - \frac{s}{M}{m_2}{\left( {N - {N_{\text{4}}}} \right)^2} \\ & \quad +\left( {\phi - \pi {m_2}} \right)\left( {N - {N_{\text{4}}}} \right)\left( {C - {C_{\text{4}}}} \right) \\ & \quad +\left( {\beta - \varepsilon } \right)\left( {C - {C_{\text{4}}}} \right)\left( {G - {G_{\text{4}}}} \right) \\ & \quad +\eta \left( {{m_1}\sigma - C} \right)\left( {C - {C_{\text{4}}}} \right)\left( {F - {F_{\text{4}}}} \right) \\ & \quad +\theta \left( {\nu {m_2} - {m_1}} \right)\left( {F - {F_{\text{4}}}} \right)\left( {N - {N_{\text{4}}}} \right)\end{aligned}$$\end{document}

choosing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${m_2}=\frac{\phi }{\pi }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${m_1}=\nu {m_2}=\frac{{\nu \phi }}{\pi }$$\end{document} we get

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}\frac{{dV}}{{dt}} & =\left[ - \frac{{p+\eta {F_4}}}{2}{{\left( {C - {C_{\text{4}}}} \right)}^2}+\eta \left( {{m_1}\sigma - C} \right) \right. \\ & \quad \left. \times \left( {F - {F_{\text{4}}}} \right)\left( {C - {C_{\text{4}}}} \right) - \frac{{{m_1}\omega }}{K}{{\left( {F - {F_{\text{4}}}} \right)}^2} \right] \\ & \quad +\left[ - \varepsilon {{\left( {G - {G_{\text{4}}}} \right)}^2}+\left( {\beta - \varepsilon } \right)\left( {C - {C_{\text{4}}}} \right) \right. \\ & \quad \left. \times \left( {G - {G_{\text{4}}}} \right) - \frac{{p+\eta {F_4}}}{2}{{\left( {C - {C_{\text{4}}}} \right)}^2} \right] - \frac{{{m_2}s}}{K}{\left( {N - {N_{\text{4}}}} \right)^2} \\ & \leq \left[ - \frac{{p+\eta {F_{\text{4}}}}}{2}{{\left( {C - {C_{\text{4}}}} \right)}^2}+\eta \left( {{m_1}\sigma - {C_{\hbox{max} }}} \right) \right. \\ & \quad \left. \times \left( {F - {F_{\text{4}}}} \right)\left( {C - {C_{\text{*}}}} \right) - \frac{{{m_1}\omega }}{K}{{\left( {F - {F_{\text{4}}}} \right)}^2} \right] \\ & \quad +\left[ - \varepsilon {{\left( {G - {G_{\text{4}}}} \right)}^2}+\left( {\beta - \varepsilon } \right)\left( {C - {C_{\text{4}}}} \right) \right. \\ & \quad \left. \times \left( {G - {G_{\text{4}}}} \right) - \frac{{p+\eta {F_{\text{4}}}}}{2}{{\left( {C - {C_{\text{4}}}} \right)}^2} \right] - \frac{{{m_2}s}}{K}{\left( {N - {N_{\text{4}}}} \right)^2} \end{aligned}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dV}}{{dt}}$$\end{document} is negative when

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\{ {\begin{array}{*{20}{c}} {\frac{{{\eta ^2}K}}{{{m_1}\omega }}{{\left( {{m_1}\sigma - {C_{\hbox{max} }}} \right)}^2}<2\left( {p+\eta {F_{\text{4}}}} \right)} \\ {\frac{{{{\left( {\beta - \varepsilon } \right)}^2}}}{\varepsilon }<2\left( {p+\eta {F_{\text{4}}}} \right)} \end{array}} \right.$$\end{document}

satisfied and the intersection of these two conditions is:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox{max} \left\{ {\frac{{\eta K}}{{{m_1}\omega }}{{\left( {{m_1}\sigma - {C_{\hbox{max} }}} \right)}^2},\frac{{{{\left( {\beta - \varepsilon } \right)}^2}}}{\varepsilon }} \right\}<2\left( {p+\eta {F_{\text{4}}}} \right)$$\end{document}

Now, we observe that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{dV}}{{dt}}$$\end{document} satisfies negative definiteness inside the attraction region ‘X’, provided that condition (38) is met.

Parameter estimation

For parameter estimation, China’s CO₂ emission data is utilized, encompassing emissions from fossil fuel combustion and land-use change while excluding other carbon emission sources [5], GDP [36], human population [37] and forest area [38] from 2000 to 2022. The natural rate of growth of atmospheric CO_2_ we take 1.68 ppm per year [39]. For the 2000–2022 period, the average annual per capita anthropogenic CO_2_ emission rate is 1.46 billion tons, equivalent to 0.08 parts per million (ppm) per million people annually. As same as [40] we take \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =0.000{\text{8}}$$\end{document} . From 2001 to 2010, the atmospheric lifetime of carbon dioxide typically ranges from 30 to 95 years [41], accordingly, the natural sink rate of atmospheric CO_2_ is 0.016. The average growth rate of GDP \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} in China from 2000 to 2022 is 0.02145. According to [38], the growth rate of forest area \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} we take 0.06133. The average intrinsic growth rate of population is approximately 0.00529 during this period of time. We take the mortality rate of the population caused by global warming as 0.00005 per ppm per year [42]. Based on parameter ranges commonly referenced in relevant studies on ecological-carbon cycle modeling, the deforestation rate and the rate of carbon dioxide sequestration by forests are set to 0.0004 and 0.0000001, respectively. Forests absorb carbon dioxide through photosynthesis to promote themselves at rate 0.01 per million hectares per year [43]. We assume that the rate of human population increase due to forest is proportional to the consumption of forest resources by a rate 0.001 [44]. As same as [40] we take \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =0.00{\text{03}}$$\end{document} . Based on the above assumptions and conjectures, the following parameter values will be used in the subsequent numerical simulation.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha ={\text{1}}{\text{.68}}$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi =0.00{\text{8}}$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =0.000{\text{3}}$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta ={\text{0}}{\text{.0000001}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=0.016$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu =0.02145$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon ={\text{0}}{\text{.0008}}$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega =0.06133$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s=0.00529$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M={\text{1720}}$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K={\text{11000}}$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta {\text{=}}0.000{\text{4}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma ={\text{0}}{\text{.01}}$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu ={\text{0}}{\text{.001}}$$\end{document} ,  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi ={\text{0}}{\text{.00005}}$$\end{document}

Numerical simulation

For the confirmation and graphical representation of analytical findings, we simulated model system (5) with the parameter values listed in Sect. 5. The numerical simulation was carried out utilizing MATLAB R2023a. The interior equilibrium components are obtained as: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_{\text{4}}}(130.9959,26.8125,{\text{3607}}{\text{.3559}},{\text{59}}{\text{.5748}})$$\end{document} . The Jacobian matrix eigenvalues at equilibrium \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_4}$$\end{document} are calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- 0.0{\text{084}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- 0.0{\text{102}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- {\text{0}}{\text{.0262}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- 0.0008$$\end{document} , all of which are negative, thus establishing the local asymptotic stability of equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_4}$$\end{document} . For the data mentioned above, the solution trajectories of the model system (5) have been plotted in Fig. 3 and Fig. 4 with different initial conditions. As observed, all trajectories starting inside the region of attraction tend toward equilibrium point \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_4}$$\end{document} , demonstrating the nonlinear stability of interior equilibrium in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C - G - F$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N - F - C$$\end{document} spaces.Fig. 3. Global Stability of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_4}$$\end{document} in C-G-F space

Fig. 4. Global Stability of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${E_4}$$\end{document} in N-F-C space

Figs. 5 and 6 depict highly significant results for the investigated dynamical system. These graphs are plotted to observe the temporal changes in the concentration of carbon dioxide \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C(t)$$\end{document} , human population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N(t)$$\end{document} and forest area \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(t)$$\end{document} . An examination of the variations is performed with regard to various parameter values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} . By comparison, the values of other parameters stay unchanged, as demonstrated in Table 1. From Fig. 5, we observe that as the anthropogenic carbon dioxide emission rate ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} ) increases from 0.005 to 0.006, the carbon dioxide concentration rises from 122.1270 ppm to 128.1255 ppm, while the forest area rises from 1569.5025 million hectares to 1615.2463 million hectares. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} increases further from 0.006 to 0.007, the CO₂ concentration increases to 134.0093 ppm, and the forest area increases to 1707.6322 million hectares. Consistent trends in CO₂ concentration and human population persist with additional increases in the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} .

From Fig. 6, as the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} increases from 0.0005 to 0.0006, the carbon dioxide concentration drops from 139.7815 ppm to 137.1956 ppm, while the human population decreases from 122.3574 million to 111.6727 million. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} further increases from 0.0006 to 0.0007, the CO₂ concentration decreases to 134.7915 ppm, and the human population drops to 93.4073 million. Consistent patterns in CO₂ concentration and human population persist with additional increments in parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} . Thus, according to these numerical findings, it can be deduced that the rate of anthropogenic carbon dioxide emissions elevates atmospheric carbon dioxide concentration and exerts a certain level of detriment to human health. Conversely, forests are capable of absorbing carbon dioxide via photosynthesis, a process that not only fosters their own growth but also markedly decreases the atmospheric carbon dioxide concentration.Fig. 5. Time series graph of C and F for different value of parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document}

Fig. 6. Time series graph of C and N for different value of parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document}

Sensitivity analysis

Global sensitivity methods such as Sobol indices and Morris screening exhibit significant advantages in detecting higher-order parameter interactions, while also holding great importance for improving sensitivity analysis [45-48]. The core focus of this study is to quantitatively analyze the interaction relationships among the four variables (atmospheric carbon dioxide, GDP, forests, and population) based on nonlinear differential equations. In this research context, PRCC is a mainstream and commonly used method for conducting parameter sensitivity analysis in nonlinear differential equation models. Its core advantage lies in its ability to accurately identify the key parameters that have a significant impact on the critical dynamic behaviors of the system. As documented by Bidah et al.[49], this methodology facilitates the evaluation of individual parameter fluctuations on the aggregate model response. A positive PRCC value denotes a direct dependency between model parameters and their outputs, such that an increment in parameter values typically elicits a pronounced rise in model output, whereas a decrement generally results in output reduction. Conversely, Fanuel et al.[50] illustrate that a negative Partial Rank Correlation Coefficient (PRCC) value indicates an inverse relationship: when the absolute value of a parameter increases, the model output decreases accordingly; conversely, when the absolute value of the parameter decreases, the model output increases.

The calculation formula of PRCC is as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & PRCC\left( {{X_i},Y} \right) \\ & \quad =\frac{{ - \sum\nolimits_{{t=1}}^{n} {\left( {{e_{Y|{X_{ - i,t}}}} \cdot {e_{X|{X_{ - i,t}}}}} \right)} }}{{\sqrt {\sum\nolimits_{{t=1}}^{n} {e_{{Y|{X_{ - i,t}}}}^{2}} } \cdot \sqrt {\sum\nolimits_{{t=1}}^{n} {e_{{Y|{X_{ - i,t}}}}^{2}} } }} \end{aligned}$$\end{document}

where: n is the total number of samples; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${X_i}$$\end{document} denotes the i-th input variable; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${X_{ - i}}$$\end{document} represents all input variables except; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e_{Y|{X_{ - i,t}}}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${e_{X|{X_{ - i,t}}}}$$\end{document} are respectively the residuals obtained by linearly regressing the rank of the output variable and the rank of the i-th input variable on the ranks of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${X_{ - i}}$$\end{document} ; PRCC ranges from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left[ { - 1,1} \right]$$\end{document} : the larger the absolute value, the stronger the variable sensitivity; positive and negative signs indicate positive and negative correlations, respectively. The baseline values and intervals of the parameters used for sensitivity analysis in this study are provided in Table 1.Fig. 7. The PRCC value of each parameter with respect to compartment C(t), G(t), F(t) and N(t) at t = 4000 respectively

When the sample size is 1000, the results of the Partial Rank Correlation Coefficient (PRCC) sensitivity analysis are illustrated in Fig.9 From Fig. 7(a), parameters exerting positive effects on the compartment \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C(t)$$\end{document} are identified as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} , while p demonstrates the most pronounced negative impact on the compartment \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C(t)$$\end{document} . Fig. 7(b) reveals that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document} respectively exert the strongest positive and negative effects on compartment \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G(t)$$\end{document} , with other parameters showing negligible influence. In Fig. 7(c), parameters contributing positively to compartment \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(t)$$\end{document} include K, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} , whereas s, p, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} exhibit the most significant negative effects; other parameters have minimal impact. Fig. 7(d) shows that p, s, K, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} positively affect compartment \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\left( t \right)$$\end{document} , while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} exert negative effects, with remaining parameters having insignificant influence.

To validate the stability of the model results, this study has supplemented a convergence analysis by generating a time-series plot of Partial Rank Correlation Coefficient (PRCC) sensitivity over a time step range from t = 1 to t = 4000 in Fig. 8. This plot was used to track the variation trend of sensitivity indices throughout the system’s operation.

As shown in the time-series plot, minor fluctuations in PRCC sensitivity were observed during the early stage of simulation, which is consistent with the dynamic adjustment characteristics of the system in its initial operational phase. With the increase in time steps, as the system gradually approached a steady state, the PRCC sensitivity indices became notably stable and no significant fluctuations were observed thereafter. This indicates that the sensitivity results of the model tend to be stable.The link between parameter estimation uncertainty and sensitivity propagation can also be elucidated through Fig. 8. For each output indicator, the fluctuation range of the PRCC values for nearly all parameters is consistently maintained during long-term evolution. This quantitative characteristic corresponds to a low impact magnitude of the estimation uncertainty of most parameters on the model outcomes through sensitivity propagation.Fig. 8. The variation of PRCC value of the parameter with respect to (a) carbon dioxide, (b) GDP, (c) forest area and (d) human population

Table 1. Parameter with their BaslineParamaterDescriptionBaslineInterval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} Natural CO₂ emission rate1.68[1.521 1.848] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document} Anthropogenic CO_2_ emission rate0.008[0.0072 0.0088] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document} GDP-driven CO₂ decay rate0.0008[0.00072 0.00088] p CO₂ decay rate0.016[0.0144 0.0176] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} CO₂ depletion coefficient due to forest area0.000001[0.0000009 0.0000011] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} GDP growth rate0.02145[0.019305 0.023595] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} Inherent growth rate of forest area0.06133[0.055197 0.067463] K Forest area carrying capacity11,000[10000 12000] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} Rate of forest loss0.0004[0.00036 0.00044] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} Growth rate of forest area for CO₂ absorption0.01[0.009 0.011] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} CO₂ emissions in the Process of GDP Growth0.0003[0.00027 0.00033] s Inherent growth rate of human population0.00529[0.004761 0.005819] M Human population carrying capacity1720[1542 1892] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} Forest area-related human population growth rate0.001[0.0009 0.0011] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi$$\end{document} Mortality rate coefficient from elevated CO₂0.00005[0.000045 0.000055]

Optimal control

The rise in carbon emissions is intrinsically linked to China’s rapid GDP growth. However, once economic development reaches a certain threshold, these economic resources can be effectively harnessed to implement technological interventions and concerted efforts aimed at absorbing and mitigating atmospheric carbon dioxide concentrations, we denote these measures by u. Nonetheless, a significant budgetary allocation is needed to fund the expenses related to these measures. Therefore, in terms of project implementation, it is necessary to formulate a cost-optimal intervention strategy, with an implementation speed that is sufficient to carry out adequate measures and actions while also minimizing the implementation cost. Given that u is not a constant but a Lebesgue measurable function over the finite time interval [0, tf​], the model system can be rewritten as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \frac{{dC}}{{dt}}=\alpha +\phi N+(\beta - u\left( t \right))G - \eta CF - pC \hfill \\ & \frac{{dG}}{{dt}}=\mu - u\left( t \right)G \hfill \\ & \frac{{dF}}{{dt}}=\omega F\left( {1 - \frac{F}{K}} \right) - \theta NF+\eta \sigma CF \hfill \\ & \frac{{dN}}{{dt}}=sN\left( {1 - \frac{N}{M}} \right)+\theta \nu NF - \pi CN \hfill \\ \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C\left( 0 \right)={C_0} \geqslant 0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G\left( 0 \right)={G_0} \geqslant 0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\left( 0 \right)={F_0} \geqslant 0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\left( 0 \right)={N_0} \geqslant 0$$\end{document} .

To minimize the objective cost function, we employ Pontryagin’s maximum principle [51]. The specific form of the objective cost function is detailed as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=\mathop {\hbox{min} }\limits_{u} \int_{0}^{{{t_f}}} {\left[ {AC\left( t \right)+\frac{B}{2}{u^2}\left( t \right)} \right]} dt$$\end{document}

Here, the parameters A and B represent the weighting parameters of the function (40). The term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{B}{2}{u^2}\left( t \right)$$\end{document} characterizes the cost incurred by the measures and actions. Subject to model (39), we seek optimal control \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u_*}(t)$$\end{document} such that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J\left( {{u_*}(t)} \right)=\mathop {\hbox{min} }\limits_{{u(t) \in \Theta }} J\left( {u(t)} \right)$$\end{document}

where the control set is denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta =$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\{ {u\left( t \right):0 \leqslant u\left( t \right) \leqslant {u_{\hbox{max} }}{\text{ for t}} \in \left[ {0,{t_f}} \right]} \right\}$$\end{document} .

Theorm 8.1

On a fixed interval [0,tf], an optimal control \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u_*} \in \Theta$$\end{document} exists to minimize the objective function (40) under the constraint of system (39).

Referring to [52], the optimal control problem under our consideration should comply with the conditions as follows:

1. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x=\left( {C(t),G(t),F(t),N(t)} \right)$$\end{document} , for a given initial value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x_0}$$\end{document} , the set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\{ {\left( {{x_0},u} \right)} \right\}$$\end{document} composed of control variable u and the solutions of the state equation that satisfy the initial conditions are non-empty.
2. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta$$\end{document} should be closed and convex. System (39) is a function of the control variable u, and the coefficients of the objective function depend on time and state variables.
3. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D=AC(t)+\frac{B}{2}{u^2}(t)$$\end{document} is convex on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta$$\end{document} and satisfies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D \geqslant f(u)$$\end{document} , here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(u)$$\end{document} is a continuous function and satisfies condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {\lim }\limits_{{\left| u \right| \to \infty }} \frac{{f(u)}}{{\left| u \right|}}=\infty$$\end{document} . Note that | · | represent the norm.

Proof According to the proof of the boundedness of the model, we know that x is bounded, that is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x_{\hbox{max} }}=\left( {{C_{\hbox{max} }},{G_{\hbox{max} }},{F_{\hbox{max} }},{N_{\hbox{max} }}} \right)$$\end{document} . As long as u is bounded within \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta$$\end{document} , the solutions of system (39) are always bounded. Therefore, condition one is satisfied.

According to the definition, given a control set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u \in \left[ {0,1} \right]$$\end{document} , thus set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta$$\end{document} is closed. According to the definition of a convex set, let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta$$\end{document} be a set. If for any \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${x_1},{x_2} \in \Theta$$\end{document} and any real number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta \in \left[ {0,1} \right]$$\end{document} , we have that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta {x_1}+(1 - \delta ){x_2} \in \Theta$$\end{document} , then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta$$\end{document} is called a convex set. Therefor, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta {x_1}+(1 - \delta ){x_2} \in \Theta$$\end{document} implying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Theta$$\end{document} is convex. System (39) can obviously be expressed as a function of the control variable u. The coefficients of the objective function A and B depend on time t and state variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C\left( t \right)$$\end{document} . Hence, we fulfilled 2.

The integrand \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D=AC(t)+\frac{B}{2}{u^2}(t)$$\end{document} is convex due to the quadratic form of u. Further more, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D=AC(t)+\frac{B}{2}{u^2}(t) \geqslant \frac{B}{2}{u^2}(t)=f(u)$$\end{document} . Obviously, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(u)$$\end{document} is continuous and satisfies condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathop {\lim }\limits_{{\left| u \right| \to \infty }} \frac{{f(u)}}{{\left| u \right|}}=\infty$$\end{document} . Hence, we fulfilled 3. Therefore, an optimal control \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u_*} \in \Theta$$\end{document} exists to minimize the objective function (40) under the constraint of system (39) over the fixed interval [0, tf].

Employing Pontryagin’s maximum principle to characterize the optimal control, the Hamiltonian is given in the following way.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} H & \left( {C,G,F,N,u,{u_1},{u_2},{u_3},{u_4}} \right)\\ & = AC(t)+\frac{B}{2}{u^2}(t) \\ & \quad +{u_1}\left[ {\alpha +\phi N+\left( {\beta - u} \right)G - \eta CF - pC} \right] \hfill \\ & \quad +{u_2}\left( {\mu - uG} \right) \hfill \\ & \quad +{u_3}\left[ {\omega F\left( {1 - \frac{F}{K}} \right) - \theta NF+\eta \sigma CF} \right] \hfill \\ & \quad +{u_4}\left[ {sN\left( {1 - \frac{N}{M}} \right)+\theta \nu NF - \pi CN} \right] \hfill \\ \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u_i}\left( {i=1,2,3,4} \right)$$\end{document} are the adjoint variable to be determined by solving the following equations.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{u^{\prime}}_1} = - \frac{{\partial H}}{{\partial C}} & = - A+{u_1}\left( {\eta F+p} \right) - {u_3}\eta \sigma F+{u_4}\pi N \hfill \\ {{u^{\prime}}_2} = - \frac{{\partial H}}{{\partial G}} & ={u_1}\left( {u - \beta } \right)+{u_2}u, \hfill \\ {{u^{\prime}}_3} = - \frac{{\partial H}}{{\partial F}} & ={u_1}\eta C - {u_3}\left[ {\omega \left( {1 - \frac{{2F}}{K}} \right) - \theta N+\eta \sigma C} \right] \\ & \quad - {u_4}\theta \nu N \hfill \\ {{u^{\prime}}_4} = - \frac{{\partial H}}{{\partial N}} & = - {u_1}\phi +{u_3}\theta F - {u_4}\\ & \quad \times \left[ {s\left( {1 - \frac{{2N}}{M}} \right)+\theta \nu F - \pi C} \right] \hfill \\ \end{aligned}$$\end{document}

along with transversality conditions

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u_1}\left( {{t_f}} \right)={u_2}\left( {{t_f}} \right)={u_3}\left( {{t_f}} \right)={u_4}\left( {{t_f}} \right)=0$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u_*}=\hbox{max} \left\{ {0,\hbox{min} \left\{ {{u_{\hbox{max} }},\frac{{\left( {{u_1}+{u_2}} \right)G}}{B}} \right\}} \right\}$$\end{document}

Numerical simulation of optimal control

To demonstrate the optimal mitigation strategies for the control of future CO_2_ level, the optimality system (39) is solved numerically by choosing the upper limit of the control \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${u_{\hbox{max} }}={\text{0}}{\text{.008}}$$\end{document} , weight parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A=0.{\text{000}}1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=10$$\end{document} , final time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${t_f}=100$$\end{document} . We utilized the forward-backward sweep method to numerically solve the optimality system corresponding to the parameter values in Sect. 5. First, we initialized the control variable with reasonable guesses. The state equations were then integrated forward in time using the fourth-order Runge-Kutta method, while the adjoint equations were solved backward in time. The control was updated iteratively until convergence. This process was repeated until the attainment of the desired convergence. Initial states are set as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C_0}=130$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G_0}=0.121$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${F_0}=1003$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${N_0}=80$$\end{document} . The solution trajectories for the concentration of carbon dioxide \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C\left( t \right)$$\end{document} , the forest area \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F(t)$$\end{document} , and the human population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N(t)$$\end{document} , both under the dynamic optimal control and in the absence of control strategies, are depicted in Fig.9. This figure clearly shows the significant reduction in carbon dioxide under the time-dependent optimal control. It is plainly evident from these graphs that the optimal control strategy outperforms the strategy without control, effectively demonstrating its superiority. The solution trajectories for atmospheric CO₂ concentration \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C\left( t \right)$$\end{document} , GDP \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G\left( t \right)$$\end{document} , forest area \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\left( t \right)$$\end{document} , and human population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\left( t \right)$$\end{document} -under both dynamic optimal control and uncontrolled conditions-are shown in Fig 0.10. Notably, the plots clearly demonstrate that the optimal control strategy outperforms the uncontrolled scenario, effectively confirming its superiority.Fig. 9. Graph trajectories with and without optimal control for: (a) carbon dioxide, (b) GDP, (c) forest area, (d) human population

Conclusion

With the rapid advancement of human society, massive greenhouse gas emissions (predominantly CO₂) have driven global warming and climate change, posing severe challenges to sustainable development. Understanding the dynamic behavior of atmospheric carbon dioxide is crucial for mitigating these environmental challenges. China currently contributes around one-third of global carbon emissions, a scale closely associated with its rapid GDP growth. As such, the study of China's carbon emission trends is of great significance for global climate change mitigation efforts. In this research, we put forward and analyze a nonlinear mathematical model that establishes a correlation between CO₂ emissions and GDP, forest area, as well as population size.

As the world’s largest developing country and a contributor to approximately one-third of global carbon emissions, China’s emission trends are closely tied to its rapid GDP growth. Given that fossil fuel combustion and industrialization (key drivers of GDP growth) account for 90% of global carbon emissions, the study of China’s GDP-CO₂ interaction is particularly significant for global climate action. Against this backdrop, this research proposed and analyzed a nonlinear mathematical model integrating four core factors: CO₂ emissions, GDP, forest area, and population size.

This model serves as a valuable tool for predicting the long-term impact of China’s GDP on atmospheric CO₂ evolution and has undergone rigorous theoretical validation: the boundedness of the system was verified using the comparison theorem; the asymptotic stability conditions of four equilibrium points were derived via Jacobian matrix eigenvalues and the Hurwitz criterion; and the global stability condition of the coexistence equilibrium point was obtained by constructing a Lyapunov function. Numerical simulations, based on the parameter values in Table 1, further confirmed the model’s practical relevance-while GDP growth is accompanied by increased carbon emissions, the model reveals that GDP can conversely be leveraged to implement mitigation measures.

To identify the key regulatory targets for such measures, this study conducted a Partial Rank Correlation Coefficient (PRCC) sensitivity analysis. The results show that the parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document} exhibit insignificant sensitivity to the system. This implies that, within the GDP-driven mitigation framework (e.g., investment in green technologies, implementation of carbon emission reduction policies), the processes represented by these two parameters have a relatively limited regulatory effect on atmospheric CO₂ concentrations. Therefore, future relevant strategies should prioritize parameters with higher sensitivity to maximize the efficiency of GDP-oriented carbon emission reduction measures.

In practice, achieving a win-win scenario of economic growth and carbon reduction requires collaborative efforts from the government, enterprises, and society at large. Targeted actions—including technological innovation, industrial structure adjustment, energy structure optimization, policy guidance, and improved public low-carbon awareness—will drive China’s transition to a low-carbon economy, ultimately contributing to both national sustainable development and global climate change mitigation.