Optimal Control of Two-Species Mutualistic Interaction under Holling Type II Functional Response

Mutualistic interactions, in which two species benefit from their biological association, are fundamental components of ecosystem function and resilience. Such interactions can be observed across a wide variety of natural systems, as shown in Figure 1, including plant-pollinator networks, mycorrhizal symbioses, and microbial cooperation. These relationships not only facilitate species coexistence but also contribute to biodiversity and ecological stability. However, the management and preservation of such mutualistic systems in the face of environmental stressors and anthropogenic disruptions remain a significant challenge in ecological science [1], [2].

In the quest to understand and predict mutualistic dynamics, mathematical models have proven to be powerful tools. Early models of mutualism often relied on simplified linear or bilinear benefit terms, which inadequately captured key ecological constraints such as resource limitations and time delays in interaction benefits. These models frequently resulted in unrealistic predictions, such as unbounded population growth [4], [5]. To address these shortcomings, we have incorporated nonlinear functional responses, particularly Holling-type, that reflect saturation effects and handling time, to provide more biologically accurate descriptions of mutualistic benefits. The Holling Type II functional response is widely used in ecological and mutualism models to introduce biological realism by incorporating saturation in interaction benefits. It prevents unbounded population growth by ensuring that the per-capita benefit decreases as population density increases, leading to a stable equilibrium instead of the unlimited growth predicted by linear Lotka-Volterra mutualism models. The parameter $h$ represents handling time, which captures biological constraints such as limited interaction time or physical capacity, and defines the maximum achievable benefit (saturation level) [6]. While this saturation stabilizes population dynamics at high densities, it can also produce Allee-type thresholds at low densities, where mutualistic interactions become insufficient to sustain growth. Because the parameters a (interaction rate) and h (handling time) can be estimated from empirical observations, the Holling Type II response provides a bridge between theoretical population models and real ecological data, making the model both analytically stable and biologically interpretable.

Several studies have advanced mutualism modeling by integrating nonlinearities and analyzing equilibrium behavior, local stability, and bifurcations. Despite the advancement, there remains a gap in the integration of mathematical rigor with practical management strategies [4], [7], [8]. In particular, few models consider implementing control strategies that balance ecological stability with resource efficiency, and even fewer examine how optimal interventions can be applied dynamically over time [9].

In this study, we developed a refined two-species mutualism model that incorporated Holling-type functional responses to reflect saturating mutualistic effects and biologically realistic interactions. The model also includes self-limiting population terms to ensure bounded growth. We first established the biological feasibility of the system by proving that its solutions are non-negative and bounded. We then performed a comprehensive qualitative analysis of the model, to examine the existence and stability of equilibrium points. This was complemented by numerical simulations under varying ecological parameters to illustrate and validate the analytical findings [10]. To assess the influence of model parameters, we conducted sensitivity analysis with a focus on how changes in key ecological factors affected critical population thresholds and system behavior [11], [12]. More importantly, our study introduced an optimal control framework to the mutualism model. We defined time-dependent control variables aimed at promoting species coexistence while minimizing ecological or economic costs. Using the fourth-order Runge–Kutta (RK4) method , we solved the optimal control problem and derived management strategies that enhanced mutualistic benefits under constrained conditions [13].

The first attempt to mathematically explore mutualistic interactions between species was done by May (1976), who introduced a system of coupled nonlinear differential equations that captures the population dynamics under mutualism. The classical mutualism model takes the following form [14]:

\[\begin{aligned} & \frac{d N_1}{d t}=\frac{r_1}{K_1} N_1\left(K_1-N_1+\alpha_{12} N_2\right) \\ & \frac{d N_2}{d t}=\frac{r_2}{K_2} N_2\left(K_2-N_2+\alpha_{21} N_1\right) \end{aligned}\]

Introducing normalized mutualistic coefficients: $m_1=\frac{\alpha_{12}}{K_1}, m_2=\frac{\alpha_{21}}{K_2}$, yields the following simplified form:

The coefficients $m_1$ and $m_2$ represent the per capita mutualistic effects: $m_1$ is the effect of one individual of

species 2 on the growth of species 1 , scaled by $K_1, m_2$ is the effect of one individual of species 1 on the growth of species 2 , scaled by $K_2$. To nondimensionalize the Eq. (1), we introduced the transformations $\widetilde{N}_1 = \frac{N_1}{K_1}$, $\widetilde{N}_2 = \frac{N_2}{K_2}$, $\tilde{t} = r_1 t$, $r = \frac{r_2}{r_1}$, $s_1 = m_1 K_2$ and $s_2 = m_2 K_1$ Thus, we obtained the non-dimensionalized system in Eq. (2):

The parameters $s_1 = m_1 K_2$ and $s_2 = m_2 K_1$ are the scaled mutualism coefficients, representing the mutualistic benefit each species provides to the other, adjusted by the respective carrying capacities. Specifically, $s_1$ measures the mutualistic benefit provided by species-2 to species-1 (accounting for species-2's carrying capacity), and $s_2$ measures the mutualistic benefit provided by species-1 to species-2 (accounting for species-1's carrying capacity) and $r = \frac{r_2}{r_1}$. From Eq. (2), it can be observed that the solutions of the system become unbounded over time for $s_1 = s_2 = 1$, as shown in Figure 2. This implies that population densities can grow without limit, leading to an outcome that is biologically unrealistic. In an ecological system, unregulated growth is constrained by various limiting factors, including resource availability, environmental carrying capacity, and intraspecific and interspecific competition. This unbounded behavior reveals a fundamental shortcoming of the current model; it fails to capture the regulatory mechanisms inherent in real ecosystems. To improve ecological realism, it is necessary to incorporate additional components that ensure bounded system dynamics, as illustrated in Figure 3.

To overcome the limitations of the classical model, we introduced a Holling Type II functional response [15] to modify the model in Eq. (2) underlying the following assumptions:

1. Two interacting mutualistic species with population densities $N_1(t), N_2(t) \geq 0$.

2. In the absence of interaction, each species follows logistic growth with intrinsic growth rates $r_1, r_2 > 0$.

3. Mutualistic interaction is described by a Holling Type II functional response with parameters $s_1, s_2, h_1, h_2 > 0$.

4. The parameter $r=\frac{r_2}{r_1}$ and since $r_1, r_2>0$, we assumed $r>0$. The case $r=0$ (i.e., $r_2=0$ ) is biologically unrealistic in the present framework.

To improve ecological realism, it is necessary to incorporate Holling Type II functional response that enforces bounded dynamics of mutualism. This adjustment improves biological realism by capturing the saturation effect in mutualistic interactions:

Here $\frac{s_1 N_2}{1+h_1 N_2}$ and $\frac{s_2 N_1}{1+h_2 N_1}$ represent saturating mutualistic benefits, thus limiting growth at high population densities, $h_1$ and $h_2$ are handling time constants that reduce mutualistic efficiency at higher densities [16]. This formulation imposes natural constraints on population growth, to ensure boundedness and eliminate biologically unrealistic outcomes.

This section is devoted to the study of the basic qualitative properties of the model, specifically the non-negativity and boundedness of its solutions.

Theorem 1. The initial value problem (IVP) as defined by the Eq. (3) with non-negative initial conditions: $N_1(0) \geq 0, N_2(0) \geq 0$. Then, the corresponding solutions $N_1(t)$ and $N_2(t)$ remain non-negative for all $t \geq 0$; that is, $N_1(t) \geq 0$ and $N_2(t) \geq 0$ for all $t \geq 0$.

Proof of Theorem 1. Assume that $N_1(0) \geq 0, N_2(0) \geq 0$, and define the following functions: \[f\left(N_1, N_2\right)=1-N_1+\frac{s_1 N_2}{1+h_1 N_2}\quad \text {and} \quad g\left(N_1, N_2\right)=1-N_2+\frac{s_2 N_1}{1+h_2 N_1}\] Then the system becomes, \[\frac{d N_1}{d t}=N_1 \cdot f\left(N_1, N_2\right) \quad \text {and} \quad \frac{d N_2}{d t}=r N_2 \cdot g\left(N_1, N_2\right)\] these yield, \[N_1(t)=N_1(0) e^{\int_0^t f\left(N_1, N_2\right) d S}\] since $N_1(0) \geq 0, N_2(0) \geq 0$, it follows that $N_1(t) \geq 0$ and $N_2(t) \geq 0$ for all $t \geq 0$. So, for any non-negative initial populations, the model guarantees non-negative solutions for all future times. Thus, the system respects the biological constraint that population densities cannot be negative. This result follows from standard methods for proving non-negativity in biological population models [17]. Q.E.D.

Theorem 2. For any initial condition $\left(N_1(0), N_2(0)\right) \in \mathbb{R}_{+}^2$, the solutions of the system (3) are uniformly bounded and remain within the region $\Omega$, where $\Omega=\left\{\left(N_1, N_2\right) \in \mathbb{R}_{+}^2 \left\lvert\, N_1 \leq 1+\frac{s_1}{h_1}+\epsilon\right., N_2 \leq 1+\frac{s_2}{h_2}+\epsilon\right\}$ for some $\epsilon > 0$ and all $t \geq t_0$.

Proof of Theorem 2. To demonstrate the solutions of the system (3) remain bounded and to determine an invariant region $\Omega$, we defined the total population function: $V=N_1+N_2$, so that, \[\frac{d V}{d t}=N_1\left(1-N_1+\frac{s_1 N_2}{1+h_1 N_2}\right)+r N_2\left(1-N_2+\frac{s_2 N_1}{1+h_2 N_1}\right).\] Holling Type II functional response is of the form $f(N) = \frac{a N}{1+b N}$ where $a > 0, b > 0$, and $N \geq 0$. This function represents a saturating response, commonly used in ecological modeling (e.g., Holling Type II functional response) [18]. As $N \rightarrow \infty$, we found $\frac{a N}{1+b N}=\frac{a}{b}$ that implied $\frac{a}{b} \geq 0$ which is essentially true, since $a, b>0$. We now estimated the upper bounds for the functional response terms $\frac{s_1 N_2}{1+h_1 N_2} \leq \frac{s_1}{h_1}, \frac{s_2 N_1}{1+h_2 N_1} \leq \frac{s_2}{h_2}$. Using these bounds, we obtained: \[\frac{d V}{d t}=N_1\left(1-N_1+\frac{s_1}{h_1}\right)+r N_2\left(1-N_2+\frac{s_2}{h_2}\right)\] Let us define $C_1=1+\frac{s_1}{h_1}$ and $C_2=1+\frac{s_2}{h_2} \cdots$ Then $\frac{d V}{d t} \leq N_1\left(C_1-N_1\right)+r N_2\left(C_2-N_2\right)$. Since $N_1\left(C_1-N_1\right)$ and $N_2\left(C_2-N_2\right)$ are concave quadratic functions with maximum values at $N_1=\frac{C_1}{2}$ and $N_2=\frac{C_2}{2}$, this implies that for large enough values of $N_1$ or $N_2$, the growth rate becomes negative. Thus, there exists a sufficiently large constant $M > 0$ such that whenever $N_1+N_2>M$, we have $\frac{d V}{d t}<0$. If $N_1+N_2$. becomes too large, then the negative growth term ensures that the total population decreases, thereby reducing $N_1+N_2$. This implies that the total population $V(t)=N_1(t)+N_2(t)$ cannot grow indefinitely and is uniformly bounded over time. Let $C = \max \left\{1+\frac{s_1}{h_1}, 1+\frac{s_2}{h_2}\right\}$ when $ N_i \geq C+\epsilon$ $\Rightarrow \frac{d V}{d t} \leq-\epsilon\left(N_1+N_2\right)<0$. For $N_1>1+\frac{s_1}{h_1}+\epsilon$, \[\frac{d N_1}{d t}<0 \quad (\text{Population decreases}).\] For $N_2>1+\frac{s_2}{h_2}+\epsilon$, \[\frac{d N_2}{d t}<0 \quad (\text{Population decreases}).\] Therefore, the solutions remain within a bounded region $\Omega \subset \mathbb{R}_{+}^2$, and the system does not exhibit unbounded growth, $\Omega=\left\{\left(N_1, N_2\right) \in \mathbb{R}_{+}^2 \left\lvert\, N_1 \leq 1+\frac{s_1}{h_1}+\epsilon\right., N_2 \leq 1+\frac{s_2}{h_2}+\epsilon\right\}$ for some $\epsilon >0$ and all $t \geq t_0$. This boundedness result follows from standard qualitative analysis techniques in population dynamics [4-19]. Q.E.D.

The four steady states of the system are: $E_1 = (0,0), E_2 = (0,1), E_3 = (1,0), E_4 = \left(N_1^*, N_2^*\right) = \left(\frac{2 a-\left(h_1+s_1\right) b \pm\left(h_1+s_1\right) \sqrt{\Delta}}{2 a-h_1 b \pm h_1 \sqrt{\Delta}}, \quad \frac{-b \pm \sqrt{\Delta}}{2 a}\right)$, where $\Delta = b^2-4 a c$ and the coefficients $a, b, c$ are as defined as

\[\begin{gathered} a=-\left(h_1+h_2\left(h_1+s_1\right)\right), \\ b=h_1+h_2\left(h_1+s_1\right)-1-h_2+s_2\left(h_1+s_1\right), \\ c=1+h_2+s_2 . \end{gathered}\]

The steady states of system (3) together with their feasibility conditions are presented in Table 1.

Lemma 1. The steady state (1,0) of system (3) is locally unstable.

Proof of Lemma 1. For the steady state (1,0), the Jacobian matrix takes the following form:

$J=\left(\begin{array}{cc}-1 & s_1 \\ 0 & r\left(1+\frac{s_2}{1+h_2}\right)\end{array}\right)$

Eigenvalues: $\lambda_1=-1, \lambda_2>0$ (since $s_2>0$ and $h_2>0)$. One eigenvalue is negative and the other positive; this steady state is a saddle point and is unstable. Q.E.D.

Lemma 2. The steady state (0,1) of system (3) is locally unstable.

Proof of Lemma 2. At (0,1), the Jacobian reduces to

$J=\left(\begin{array}{cc}1+\frac{s_1}{1+h_1} & 0 \\ r s_2 & -r\end{array}\right)$

Eigenvalues: $\lambda_1>0 \cdot\left(\right.$ since $s_1>0 \cdot$ and $\left.h_1>0\right), \lambda_2=-r<0 \cdot($ since $r>0)$. One eigenvalue is positive and the other negative; this steady state is a saddle point and is unstable. Q.E.D.

The non-trivial steady state $\left(N_1^*, N_2^*\right)$ of the system (3) is locally asymptotically stable if $P<0$ and $Q>0$.

undefined

When $h_1=h_2=0$, Figure 2 illustrates that this condition could result in biologically unrealistic outcomes such as unbounded population growth, as mutualistic benefits are no longer limited. To ensure ecological realism and system stability, we introduced the handling time parameters $h_1$ and $h_2$ [20]. Figure 3 displays the resulting phase portrait under these regulated mutualistic interactions. In both figures, species $N_1$ is represented by the magenta curve, and species $N_2$ by the orange curve.

Incorporating handling time produced ecologically realistic scenario [20-21], as shown in Figure 3. The solution of the system reached a stable interior equilibrium at approximately $\left(N_1^*, N_2^*\right) \approx(10.1,10.1)$, where both species coexisted at high densities due to mutualistic benefits. All trajectories with $N_1, N_2>0$ converged to such an equilibrium. The absence of limit cycles or extinction events reflected strong and regulated mutualism. The saturation effects introduced by $h_1$ and $h_2$ capped mutualistic gains, thus preventing unbounded growth and ensuring realistic coexistence [21].

The sensitivity index quantifies how the equilibrium population $N_i^*$ responds to changes in a parameter $p_i$ which is defined [22-24] as:

\[\gamma_{p_j}^{N_i^*}=\left(\frac{\partial N_i^*}{\partial p_j}\right)\left(\frac{p_j}{N_i^*}\right)\]

where, $p_j \in\left\{s_1, s_2, h_1, h_2, r\right\}$ and the parameters are chosen $s_1=0.5, s_2=0.5, h_1=0.3, h_2=0.2, r=1$. Table 2 summarizes the local sensitivity indices corresponding to the model parameters.

Bar diagrams below are provided to visually summarize the sensitivity of equilibrium populations $N_1^*$ and $N_2^*$ to changes in the key parameters [22] in Figure 4 and Figure 5.

To compute a numerical solution of the system (3), we considered the following IVP. We used the Runge-Kutta 4th order method in solving the IVP (4):

Figure 6, Figure 7, and Figure 8 illustrate the time series and phase portraits of two interacting species $N_1(t)$ and $N_2(t)$ under different parameter settings of $s_1, s_2, h_1, h_2$ and $r$. We observed that varying the mutualistic strength ($s_1, s_2$) significantly affected the stability and equilibrium levels of the populations. The phase portraits further confirm how initial conditions evolve toward equilibrium or diverge depending on the parameter combinations.

This section introduces a controlled version of the mutualistic system in order to investigate how external interventions can influence population persistence and ecological stability. We incorporated time-dependent control functions $q_1 (t)$ and $q_2 (t)$, representing management strategies such as species augmentation or regulated harvesting. Necessary optimality conditions are derived using Pontryagin’s Maximum Principle [25]. The controlled mutualism model is described by:

where, $q_1$ and $q_2$ represent external controls on the species populations $N_1$ and $N_2$, respectively. The ecological interpretations associated with the steady states of system (5) are listed in Table 3.

Admissible Control Set:

To ensure mathematical consistency, we defined the admissible control set as

\[U=\left\{\left(q_1, q_2\right) \in L^{\infty}(0, T): q_{i, \min } \leq q_i(t) \leq q_{i, \max }, i=1,2\right\}.\]

This formulation allows both enhancing controls $\left(q_i>0\right)$ and suppressing controls $\left(q_i<0\right)$, depending on the biological context. However, in the present instability and sustainability analysis, we focused on the biologically relevant case $q_1(t)>0, q_2(t)>0$, which corresponds to conservation or augmentation strategies.

Model Assumptions:

The model was developed under the following assumptions:

The mutualistic interaction follows a Holling Type II functional response with parameters $s_1, s_2, h_1, h_2 > 0$.,The parameter $r=\frac{r_2}{r_1}$ represents the ratio of intrinsic growth rates. Since $r_1, r_2>0$, we assumed $r>0$. The case $r=0$ (i.e., $r_2=0$ ) would eliminate the intrinsic growth of the second species and was therefore excluded.,Initial conditions satisfy $N_1(0), N_2(0) > 0$.,When we will consider $q_1(t)=q_2(t)=0$, Eq. (5) reduces to the uncontrolled mutualism model in Eq. (3).

The trivial steady state (0,0) of system (5) is locally unstable.

undefined

The semi-trivial steady state $\left(1+q_1, 0\right)$ of system (5) is locally unstable.

undefined

The semi-trivial steady state $\left(0,1+\frac{q_2}{r}\right)$ of system (5) is locally unstable.

undefined

The non-trivial steady state $\left(N_1^*, N_2^*\right)$ of system (5) is asymptotically stable if $P<0$ and $Q>0, P$ and $Q$ are defined below.

undefined

Formulating an optimal control problem in ecological systems, such as the controlled mutualism model, allows us to identify time-dependent strategies that maximize ecological benefits while minimizing associated costs. In real-world applications, resources for intervention such as conservation funding, habitat management, or species support are limited. Optimal control theory helps balance the trade-off between investing in control measures and achieving desired ecological outcomes, such as sustaining mutualistic species at healthy population levels. By designing dynamic and time-varying control policies, this approach provides a scientifically grounded method to determine when, how much, and for how long interventions should be applied [26-29]. We considered the controlled mutualism model (5) where $N_1(0) \geq 0, N_2(0) \geq 0$.

The objective of our study is to find an optimal control pair ($q_1^*$,$q_2^*$) that maximizes population levels while minimizing the costs associated with the controls over a finite time horizon $[ 0, T]$. The objective functional $J$ is defined as

\[J\left(q_1, q_2\right)=\int_0^T\left(A_1 N_1(t)+A_2 N_2(t)-\frac{1}{2} B_1 q_1^2(t)-\frac{1}{2} B_2 q_2^2(t)\right) d t\]

where, $A_1, A_2>0$ are weight constants and $B_1, B_2>0$ represent the cost of implementing the controls, $q_1$ and $q_2$ are bounded measurable control functions representing external interventions (e.g., harvesting, treatment, or stimulation).

To derive the necessary conditions for the optimal control, we applied the PMP. We defined the Hamiltonian function $H$ for the system as follows:

\[H=\left(A_1 N_1+A_2 N_2-\frac{1}{2} B_1 q_1^2-\frac{1}{2} B_2 q_2^2\right)+\lambda_1 \frac{d N_1}{d t}+\lambda_2 \frac{d N_2}{d t}\]

Substituting the expressions for the derivatives from Eq. (5)

\[\begin{aligned} H=\left(A_1 N_1+A_2 N_2\right. & \left.-\frac{1}{2} B_1 q_1^2-\frac{1}{2} B_2 q_2^2\right) \\ & +\lambda_1\left[N_1\left(1-N_1+\frac{s_1 N_2}{1+h_1 N_2}\right)+q_1 N_1+\lambda_2\left[r N_2\left(1-N_2+\frac{s_2 N_1}{1+h_2 N_1}\right)+q_2 N_2\right]\right] \end{aligned}\]

where, $\lambda_1(t)$ and $\lambda_2(t)$ are the adjoint variables. According to the PMP, the adjoint equations are given by: $\dot{\lambda}_1=-\frac{\partial H}{\partial N_1}, \dot{\lambda}_2=-\frac{\partial H}{\partial N_2} \cdot$ Differentiating $H$ with respect to the state variables $N_1$ and $N_2$, we obtained

with the transversality conditions at the final time $T:\lambda_1(T)=0$ and $T:\lambda_2(T)=0$ [28].

The optimal control pair $\left(q_1^*, q_2^*\right)$ must satisfy the optimality condition $\frac{\partial^H}{\partial q_i}=0$. Taking the partial derivatives of the Hamiltonian with respect to $q_1$ and $q_2:$

\[\begin{aligned} & \frac{\partial H}{\partial q_1}=-B_1 q_1+\lambda_1 N_1=0 \Rightarrow q_1=\frac{\lambda_1 N_1}{B_1} \\ & \frac{\partial H}{\partial q_2}=-B_2 q_2+\lambda_2 N_2=0 \Rightarrow q_2=\frac{\lambda_2 N_2}{B_2} \end{aligned}\]

Considering the admissible control set $U=\left\{\left(q_1, q_2\right): a_i \leq q_i \leq b_i\right\}$, the characterization of the optimal controls is given by:

Thus, the optimal control problem consists of solving the state system (5), the adjoint system (6), and the control characterization Eq. (7) simultaneously. This can be done numerically using the forward-backward sweep method. Solve the state system (5) forward in time from given initial conditions, starting with an initial guess for the controls $q_1(t), q_2(t)$. Using the state trajectories, solve the adjoint system (6) backward in time from the terminal conditions $\lambda_1(T)=0, \lambda_2(T)=0$. Update the control functions $q_1(t), q_2(t)$ using the control characterization (7). Repeat these steps iteratively until the controls converge within a prescribed tolerance.

The numerical results are obtained using the parameters listed in Table 4. These values are chosen to represent a mutualistic system where both species benefit from each other with a saturation effect.

Figure 9 and Figure 10 illustrate the temporal evolution of the two mutualistic populations. It is observed that the species reach a stable steady state governed by their intrinsic growth rates and mutualistic interactions. Moreover, it is observed that the incorporation of the control functions $q_1(t)$ and $q_2(t)$ significantly increases the population densities compared with their natural equilibrium levels. Consequently, optimal control strategies can effectively enhance the persistence and growth of the mutualistic species.

Table 5 defines the biological and control parameters used in our simulations. We categorized our analysis into three scenarios based on cooperation strength $s_i$, the ratio of growth rates $r$, and other parameters.

Figure 11, Figure 12, and Figure 13 illustrate the temporal dynamics of the two mutualistic populations $N_1$, $N_2$ under weak, moderate, and strong mutualistic interactions, respectively. For the strong mutualism case, both populations grew rapidly due to the high cooperative benefit and eventually reached a stable equilibrium. The application of optimal control further enhanced the growth rate and drove the system to higher equilibrium population densities compared with the uncontrolled scenario. In the case of moderate mutualism, the populations also increased and stabilized at intermediate density levels. The controlled system again showed improved growth performance, with both species attaining higher equilibrium densities than those observed in the natural system. In weak mutualism, the natural system exhibited slower growth and lower equilibrium population levels due to limited cooperative interactions between the species. However, when optimal control was applied, the populations increased more efficiently and reached noticeably higher densities. Overall, the results clearly highlighted that optimal control strategies significantly enhanced population growth and persistence across all levels of mutualistic interaction. The improvement was particularly important in the weak mutualism scenario, where external intervention was crucial for sustaining higher population densities.

This research presented a thorough exploration of optimal control strategies for managing mutualistic interactions between two ecologically interconnected species. By developing a biologically realistic mathematical model, we effectively captured the essence of mutualism, incorporating nonlinear functional responses that reflect interactions observed in natural ecosystems. We established the model’s biological feasibility by proving that the solutions are non-negative and bounded. A detailed qualitative analysis provided insight into the dynamic behavior of the system, particularly the conditions for the existence and local stability of equilibrium points. These theoretical results were reinforced by numerical simulations across a range of ecological scenarios, thus confirming the robustness and validity of the model under varying parameter values. A sensitivity analysis was conducted to assess how changes in specific parameters influenced steady-state population levels. This analysis identified the most influential ecological factors, offering valuable insights into the responsiveness and resilience of the system. To enhance the practical relevance of our model, we formulated and solved an optimal control problem by incorporating time-dependent control variables into the mutualistic framework. The results of optimal control provided guidance for sustainable population management and efficient resource allocation in mutualistic ecosystems. This study demonstrated that optimal control theory provided a powerful framework for enhancing mutualistic population dynamics. The comparative analysis revealed that properly chosen control strategies could enhance mutualistic interactions and support greater long-term population abundance. Furthermore, when control was applied, the beneficial relationship between the two species strengthened, and their populations increased. As global ecological challenges intensify, mathematically grounded management strategies become increasingly essential for maintaining biodiversity and ecosystem function. Optimal control theory plays an effective role in enhancing population growth in cooperative species systems.