*Article* **Analytical Solutions of Microplastic Particles Dispersion Using a Lotka–Volterra Predator–Prey Model with Time-Varying Intraspecies Coefficients**

**Lindomar Soares Dos Santos 1, José Renato Alcarás 1, Lucas Murilo Da Costa 1, Mateus Mendonça Ramos Simões <sup>2</sup> and Alexandre Souto Martinez 1,3,\***


**Abstract:** Discarded plastic is subjected to weather effects from different ecosystems and becomes microplastic particles. Due to their small size, they have spread across the planet. Their presence in living organisms can have several harmful consequences, such as altering the interaction between prey and predator. Huang et al. successfully modeled this system presenting numerical results of ecological relevance. Here, we have rewritten their equations and solved a set of them analytically, confirming that microplastic particles accumulate faster in predators than in prey and calculating the time values from which it happens. Using these analytical solutions, we have retrieved the Lotka–Volterra predator–prey model with time-varying intraspecific coefficients, allowing us to interpret ecological quantities referring to microplastics dispersion. After validating our equations, we solved analytically particular situations of ecological interest, characterized by extreme effects on predatory performance, and proposed a second-order differential equation as a possible next step to address this model. Our results open space for further refinement in the study of predator–prey models under the effects of microplastic particles, either exploring the second-order equation that we propose or modify the Huang et al. model to reduce the number of parameters, embedding in the time-varying intraspecies coefficients all the adverse effects caused by microplastic particles.

**Keywords:** microplastics; Lotka–Volterra; predator; prey; predator–prey; pollution; bioaccumulation; biomagnification; predation; two-species model

#### **1. Introduction**

The discovery of Bakelite in 1907 revolutionized modern life by introducing plastic materials to the world. The popularization of these polymers began with their commercial production around 1950 [1,2]. Plastic's versatility, stability, low weight, and low production cost leveraged its global market for this material [2]. The increase in plastic consumption has had serious repercussions on nature. It is estimated that 9.5 million tons of plastic end up in the oceans every year [3,4] and there are still no estimates for the amount of plastic deposited on land [5].

All discarded plastic is subject to weather effects of different ecosystems. This material can be slowly degraded by photo-oxidation, thermal pathways, mechanochemical interactions or biodegradation [6]. The results of all these processes are small particles known as microplastics, which have dimensions that vary between 1 μm and 5 mm [7] and the most different compositions, colors, and shapes [8].

**Citation:** Dos Santos, L.S.; Alcarás, J.R.; Da Costa, L.M.; Simões, M.M.R.; Martinez, A.S. Analytical Solutions of Microplastic Particles Dispersion Using a Lotka–Volterra Predator–Prey Model with Time-Varying Intraspecies Coefficients. *Math. Comput. Appl.* **2022**, *27*, 66. https://doi.org/ 10.3390/mca27040066

Academic Editor: Gabriel B. Mindlin

Received: 16 June 2022 Accepted: 1 August 2022 Published: 3 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Due to their small size, microplastics have spread across the planet. They have already been detected in all parts of the ocean [9], in the poles [10], in drinking water [11], in arable areas and pastures [12], in the habitat of terrestrial animals [13–15], and in various foods consumed by people [16]. Many studies suggest that the near-ubiquity of microplastics means that their transfer is not limited to food chain related dynamics, but also occurs through bioaccumulation from one trophic level to another [17–19].

The presence of microplastics in living organisms can have several consequences. Studies indicate the possibility of gastrointestinal inflammation [20], decreased reproductive capacity [21], reduced ability to feed [22], reduced growth rate [23], and even malformation of embryos [24]. Early research suggests that microplastics can alter the interaction between prey and predator [25], producing interesting study approaches that use Lotka–Volterra prey–predator models modified to take into account the microplastics [26].

The Lotka–Volterra model (LVM) was originally developed from the logistic equation, describing chemical reactions [27]. However, over the years, modifications of this model have allowed new and specific applications, such as the study of the interaction between prey and predators [28–30], the influence of harmful elements on population dynamics [31], and, more specifically, the effect that microplastics has on trophic relationships [18,19,26].

This paper aims to go beyond the numerical solutions of a prey–predator model by studying analytical solutions of the model proposed by Huang et al. [26], which takes into account the impact of microplastics on the population dynamics. We rewrite the Huang et al. four equations model, reducing it to two equations equivalent to one and including the timevarying intraspecies coefficients. This approach clarifies the model parameters meanings and allows us to solve analytically three special cases of ecological relationships. These special cases are based on possible extreme effects of predatory performance reduction caused by exposure to MP particles and are mathematically characterized by the decoupling of the differential equations of the model, for which we also perform numerical simulations. We also propose a second-order differential equation as a possible next step to address this model.

We organized the presentation as follows: In Section 2, we introduce Huang et al.'s work, including their model and main results. In addition, we justify our interest in their research, presenting our motivations and intentions, and rewrite four original equations of the predator–prey model, reducing them to only two ones, redefining/regrouping some quantities, including time-varying intraspecies coefficients, giving them ecological meanings analogous to the standard LVM. In Section 3, we analytically validate that microplastic (MP) particles accumulate faster in predators than in preys and calculate the characteristics times from which their concentration and changing rate of the total amount are greater in predators than in preys. After validating our model, we explore analytically and numerically three special ecological regimes characterized by extreme effects on predatory performance, which can lead these two populations to become independent. In addition, we introduce a second-order differential equation as a possible future study of our two equations system. In Section 4, we compile our results and interpretations, presenting their implications and research possibilities of refining LVM under effects of MP particles, exploring our second-order equation or modifying the standard model to reduce its number of parameters, embedding in the time varying intraspecies coefficients all the adverse effects caused by MP particles. In Section 5, we summarize our main results and point to future studies.

#### **2. Materials and Methods**

In 2020, Huang et al. [26] published a study using the LVM where they theoretically investigated predator–prey population dynamics in terms of toxicological response intensity to microplastic (MP) parts and examined the negative effects on prey feeding ability and predator performance due to MP particles. The study suggests that dynamic LVMs can be an important tool to predict the ecological impacts of MP particles on predator–prey population dynamics.

Combining a single-species model with the LVM, Huang et al. obtained the following model:

$$\begin{cases}
\dot{\mathbf{x}}\_{1}(t) = \mathbf{x}\_{1}(t)[r\_{10} - d\_{1} - r\_{11}\mathbf{C}\_{1}(t) - (a\_{1} - d\_{3})\mathbf{x}\_{2}(t)] \\
\dot{\mathbf{x}}\_{2}(t) = \mathbf{x}\_{2}(t)[-r\_{20} - r\_{21}\mathbf{C}\_{2}(t) + (a\_{2} - d\_{2})\mathbf{x}\_{1}(t)] \\
\dot{\mathbf{C}}\_{1}(t) = \mathbf{S}\_{1}\mathbf{C}\_{E} - \mathbf{g}\_{1} \\
\dot{\mathbf{C}}\_{2}(t) = \mathbf{S}\_{2}\mathbf{C}\_{E} + k\mathbf{C}\_{1}(t) - \mathbf{g}\_{2}
\end{cases} \tag{1}$$

where *x*1(*t*), *x*2(*t*), *C*1(*t*), *C*2(*t*) ≥ 0. This model neglects the influence of the intraspecific competition and shows the effect of the microplastics by its concentration. Considering *t* ≥ 0 in those differential equations, *x*<sup>1</sup> and *x*<sup>2</sup> represent the population of prey and predator, respectively. *r*10*x*<sup>1</sup> and *r*20*x*<sup>2</sup> are the intrinsic growth rate of prey and mortality rate of predator without toxicity, respectively. *a*1*x*1*x*<sup>2</sup> is the lost amount of prey eaten by predators, and *a*2*x*1*x*<sup>2</sup> is the increasing number of predators due to the feeding of prey, with *a*1, *a*<sup>2</sup> > 0. The parameters *d*1, *d*2, and *d*<sup>3</sup> represent the decline in the prey feeding ability, the adverse effect of reduced predatory performance, and the lost amount of prey eaten by the predator, respectively. The response intensities of MP particles on prey and predator are denoted by *r*<sup>11</sup> and *r*21, respectively. The microplastics egestion rates of prey and predator, *g*1, *g*<sup>2</sup> ≥ 0 (*g*<sup>1</sup> < 0 or *g*<sup>2</sup> < 0 would imply a MP particles "negative egestion", by prey or predator), are independent of the microplastics concentration in the environment *CE* ≥ 0, and the amount of microplastics concentration removed at each time step is independent of the total amount of microplastics in the organisms, *C*<sup>1</sup> and *C*2. The quantity *kC*<sup>1</sup> represents the accumulated toxicity of MP particles transferred from the prey. Finally, *S*1, *S*<sup>2</sup> ≥ 0 are related to the effects of plastic particle selection of prey and predator, respectively (*S*<sup>1</sup> < 0 or *S*<sup>2</sup> < 0, which would imply an environment MP particle "negative ingestion" rate by prey or predator).

The authors estimated parameters and performed simulations that indicated that predators are more vulnerable than prey under exposure to microplastics. The effect of MP particles on both population growths can be negligible when toxicological response intensities of prey and predator are small, the system is prey-dependent for predator functional response, and the reduced feeding capacity of prey and predator induced by microplastics does not significantly affect the population dynamics. The conclusions are compatible with empirical evidence. This study indicates that this model is adequate to approach the prediction of population dynamics of the predator–prey system under toxicological effects of persistent organic pollutants.

Stimulated by the relevance of the study and the innovativeness of its approach, we decided to explore the Huang et al. model, going beyond its numerical solutions by obtaining analytical solutions that reproduce some of their results, using these solutions to validate our equations, and going even further presenting and analyzing particular cases of possible ecological interest.

After elaborating the model, Huang et al. estimated its parameters and performed numerical simulations, implementing them using MATLAB programming and its Simulink toolbox, based on differential equations such as those presented by the system (1). We have realized that it is possible to rewrite only two equations, reducing them to the derivatives *x*˙1(*t*) and *x*˙2(*t*). In addition, it is interesting to redefine/regroup some quantities of the model, significantly reducing its number of parameters and clarifying its understanding and ecological interpretation. Since some of the parameters in the system are related, we can regroup them into effective (net) and variable rates of decline and growth.

Initially, we define the quantities: *α* = *a*<sup>1</sup> − *d*3, which is the difference between the rate of prey population decline and the rate of the lost amount of prey eaten by predators (that we call effective rate of prey population decline), and *α* = *d*<sup>2</sup> − *a*2, which is the difference between the rate of reduced predatory performance and the rate of predator population growth (that we call negative of the effective rate of predator population growth), and rewrite the system (1):

$$\begin{cases}
\dot{\mathbf{x}}\_1(t) = \mathbf{x}\_1(t)[r\_{10} - d\_1 - r\_{11}\mathbf{C}\_1(t) - a\mathbf{x}\_2(t)] \\
\dot{\mathbf{x}}\_2(t) = \mathbf{x}\_2(t)[-r\_{20} - r\_{21}\mathbf{C}\_2(t) - a'\mathbf{x}\_1(t)] \\
\dot{\mathbf{C}}\_1(t) = \mathbf{S}\_1\mathbf{C}\_E - \mathbf{g}\_1 \\
\dot{\mathbf{C}}\_2(t) = \mathbf{S}\_2\mathbf{C}\_E + k\mathbf{C}\_1(t) - \mathbf{g}\_2
\end{cases},\tag{2}$$

where *x*1(*t*), *x*2(*t*), *C*1(*t*), *C*2(*t*) ≥ 0.

We see that the solution to *C*1(*t*) is a first degree polynomial:

*C*1(*t*) = *C*1(0) + *C*˜ <sup>1</sup>*t* , (3)

with

$$
\vec{\mathbb{C}}\_1 = \mathbb{S}\_1 \mathbb{C}\_E - \mathbb{g}\_{1\prime} \tag{4}
$$

which means that the total amount of microplastics in the prey population varies at a constant rate *C*˜ 1. Note that the parameters *S*1, *CE*, and *g*<sup>1</sup> combine to form a quantity with units of [*C*1]/[*t*].

Using the solution of *C*1(*t*), we obtain the solution of *C*2(*t*), which is a second degree polynomial:

$$\mathcal{C}\_2(t) = \mathcal{C}\_2(0) + \tilde{\mathcal{C}}\_2 t + \frac{k\tilde{\mathcal{C}}\_1}{2}t^2 \,. \tag{5}$$

with

$$\mathbf{C}\_2 = \mathbf{S}\_2 \mathbf{C}\_E - \mathbf{g}\_2 + k \mathbf{C}\_1(0),\tag{6}$$

meaning that the total amount of microplastics in the predator population, unlike the prey population, varies at a rate that is directly proportional to time, and equal to *C*˜ <sup>2</sup> + *kC*˜ 1*t*.

Comparing the expressions for the total amount of MP particles in the prey, *C*1(*t*), and predators, *C*2(*t*), we show that, over time, the concentration of microplastics inside the predator population will eventually be larger than the concentration inside prey population, which occurs at

$$t\_{C\_2>C\_1} > \frac{\left(\tilde{C}\_1 - \tilde{C}\_2\right) + \sqrt{\left(\tilde{C}\_1 - \tilde{C}\_2\right)^2 - 2k\tilde{C}\_1\left(C\_2(0) - C\_1(0)\right)}}{k\tilde{C}\_1} \,. \tag{7}$$

This result analytically confirms the already known fact [26] that MP particles tend to accumulate faster in predators than in prey, which occurs when the rate of change of the total amount of microplastics inside predators is greater than the rate of change of the total amount of microplastics inside prey, more specifically at

$$t\_{\mathcal{C}\_2 > \mathcal{C}\_1} > \frac{\tilde{\mathcal{C}}\_1 - \tilde{\mathcal{C}}\_2}{k \mathcal{C}\_1} \,. \tag{8}$$

Replacing *C*1(*t*) and *C*2(*t*) in the first two equations of the system (2), one needs only two coupled differential equations to write the Huang et al. model:

$$\begin{cases}
\dot{\mathbf{x}}\_{1}(t) = \mathbf{x}\_{1}(t) \{ r\_{10} - d\_{1} - r\_{11} [\mathbf{C}\_{1}(0) + \tilde{\mathbf{C}}\_{1} t] - a \mathbf{x}\_{2}(t) \} \\
\dot{\mathbf{x}}\_{2}(t) = \mathbf{x}\_{2}(t) \left\{ -r\_{20} - r\_{21} \left[ \mathbf{C}\_{2}(0) + \mathbf{C}\_{2} t + \frac{\mathbf{k} \tilde{\mathbf{C}}\_{1}}{2} t^{2} \right] - a' \mathbf{x}\_{1}(t) \right\} \\
\end{cases} \tag{9}$$

Regrouping further this system of equations, one sees that it is possible to rearrange some of its quantities, to obtain variable rates of growth/decline of prey/predators. Regrouping the quantities: *c*<sup>0</sup> = *r*<sup>10</sup> − *d*<sup>1</sup> − *r*11*C*1(0), *c* <sup>0</sup> <sup>=</sup> <sup>−</sup>*r*<sup>20</sup> <sup>−</sup> *<sup>r</sup>*21*C*2(0), *<sup>c</sup>*<sup>1</sup> <sup>=</sup> *<sup>r</sup>*11*C*˜ 1, *c* <sup>1</sup> <sup>=</sup> *<sup>r</sup>*21*C*˜ 2, and *c* <sup>2</sup> = *<sup>r</sup>*21*kC*˜ 1/2, we define *β*1(*t*) = *c*<sup>0</sup> − *c*1*t*, which is the rate of prey population growth, and *β*2(*t*) = *c* <sup>0</sup> − *c* <sup>1</sup>*t* − *c* 2*t* 2, which is the rate of predator population decline. Note that the effect of the initial conditions appears only on the coefficients *c*<sup>0</sup> and *c* 0.

In this way, for *x*˙*i*(*t*), we have:

$$\begin{cases}
\dot{\mathbf{x}}\_1(t) = \mathbf{x}\_1(t)[\beta\_1(t) - a\mathbf{x}\_2(t)] \\
\dot{\mathbf{x}}\_2(t) = \mathbf{x}\_2(t)[\beta\_2(t) - a'\mathbf{x}\_1(t)]
\end{cases}.\tag{10}$$

This system is similar to the well-known standard Lotka–Volterra (predator–prey) equations, except that its prey population growth and predator population decline rates are now time-varying intraspecific coefficients, which have built-in the effects of microplastics in the organisms.

A simple steady-state analysis *x*˙1(*t*) = *x*˙1(*t*) = 0, as *t* → ∞, leads to the observed regimes. In one hand, if *x*1(∞) = 0, then *x*2(∞) ≥ 0 may have arbitrary values. On the other hand, if *x*2(∞) = 0, then *x*1(∞) ≥ 0 may have arbitrary values. In case *x*1(∞) = 0, then *x*2(∞) = *β*2(∞)/*α* → ∞. Similarly, *x*2(∞) = 0, then *x*1(∞) = *β*1(∞)/*α* → ∞. Particular cases presented below show these asymptotic values.

Huang et al. investigate the toxicological effects of MP particles according to the value of response intensity *r*<sup>11</sup> and *r*<sup>12</sup> implementing the model (1) using MATLAB programming and its Simulink toolbox. The authors classified interactions into four different conditions (response intensities of prey and predator to toxicological effects induced by MP particles): (a) without the influence of MP particles (*C*<sup>1</sup> = *C*<sup>2</sup> = *r*<sup>11</sup> = *r*<sup>21</sup> = 0); (b) predator and prey have the same response strength to MP particles (Δ = *r*11/*r*<sup>21</sup> = 1.0); (c) predator has much larger response strength than prey (Δ = *r*11/*r*<sup>21</sup> = 0.1); (d) predator has much smaller response strength than prey (Δ = *r*11/*r*<sup>21</sup> = 10.0). Huang et al. constructed phaseportraits and short-term population dynamics graphs of predator–prey for each of these four conditions, with the following values: *CE* = 30, *x*1(0) = 100, *x*2(0) = 100, *C*1(0) = 0, *C*2(0) = 0, *r*<sup>10</sup> = 4.1, *r*<sup>20</sup> = 4.0, *d*<sup>1</sup> = 0.1, *d*<sup>2</sup> = 0.002, *d*<sup>3</sup> = 0.002, *a*<sup>1</sup> = 0.052, *a*<sup>2</sup> = 0.052, *g*<sup>1</sup> = 1.2, *g*<sup>2</sup> = 1.3, *S*<sup>1</sup> = 0.042, *S*<sup>2</sup> = 0.039, *a*<sup>1</sup> = 0.052, *k* = 2.0. In addition, they analyzed the negative effects of PM particles on prey feeding ability and predatory performance increasing the values of *d*1, *d*2, and *d*<sup>3</sup> to *d*<sup>1</sup> = 0.6, *d*<sup>2</sup> = 0.012, and *d*<sup>3</sup> = 0.012, and taking three conditions of different response strength (*r*<sup>11</sup> = *r*<sup>21</sup> = 10.0; *r*<sup>11</sup> = 1.0; and *r*<sup>21</sup> = 10.0; and *r*<sup>11</sup> = 10.0 and *r*<sup>21</sup> = 1.0).

#### **3. Results**

Implementing the model of Equation (10) using the SciPy Python library, with the function *odeint* from *Scipy.integrate* package, and considering the five scenarios described in the previous section, we successfully reproduced the phase-portraits and the population dynamics graphs of Huang et al. The system we built consists of only two first-order coupled equations, has a reduced and more comprehensible set of parameters, and reproduces all the results of Huang et al. These equations are decoupled in three situations: *α* = *α* = 0, *α* = 0, and *α* = 0, which address specific ecological regimes [32].

#### *3.1. Maximum Reduction of Predatory Performance: α* = 0 *and α* = 0

Exposure to MP particles can cause anomalous behaviors that impair the prey feeding ability and predatory performance of organisms, which lead, for example, to a decrease in the intrinsic growth rate of prey and a reduction in the number of predators. The Huang et al. model considers some of these harmful effects. For example, it denotes *d*<sup>2</sup> to the adverse effect of reduced predatory performance, and *d*<sup>3</sup> to the consequent lost amount of prey eaten by the predator. In our study, we initially considered the extreme situation in which: 1—the presence of MP particles impairs the predatory performance to the point where there is no further increase in the predators population due to their prey feeding or, in other words, *d*<sup>2</sup> is large enough to reduce to zero the effective rate of predator population growth (−*α* = *a*<sup>2</sup> − *d*<sup>2</sup> = 0); and 2—as a consequence of the drastic reduction in the predatory performance, no more prey is eaten by predators; in other words, *d*<sup>3</sup> is large enough to reduce to zero the effective rate of prey population decline (*α* = *a*<sup>1</sup> − *d*<sup>3</sup> = 0). This is a very special scenario, where the harm caused by microplastics on predators makes these two

populations independent, transforming a relationship of predation into a relationship of neutralism, in a process described by the decoupling of the equations for *x*˙ <sup>1</sup> and *x*˙ 2.

Assuming *α* = *α* = 0, the system (10) becomes:

$$\begin{cases}
\dot{x}\_1(t) = x\_1(t)\beta\_1(t) \\
\dot{x}\_2(t) = x\_2(t)\beta\_2(t)
\end{cases} \tag{11}$$

which describes the population dynamics of two independent groups, prey, and predators, i.e., the variation in the number of individuals in one group does not affect the growth dynamics of the other.

Solving the system above, we get independent solutions for *x*<sup>1</sup> and *x*2, which will be labeled (*ind*) for future reference:

$$\mathbf{x}\_1^{(ind)}(t) = \mathbf{x}\_1(t) = \mathbf{x}\_1(0)e^{\int\_0^t dt' \beta\_1(t')} = \mathbf{x}\_1(0)e^{c\_0t - c\_1t^2/2} \tag{12}$$

and

$$\mathbf{x}\_{2}^{(ind)}(t) = \mathbf{x}\_{2}(t) = \mathbf{x}\_{2}(0)e^{\int\_{0}^{t} dt' \beta\_{2}(t')} = \mathbf{x}\_{2}(0)e^{c\_{0}'t - c\_{1}'t^{2}/2 + c\_{2}'t^{3}/3}.\tag{13}$$

Analyzing the Equations (12) and (13), it is convenient to calculate *X*(*ind*) <sup>1</sup> (*t*) and *X*(*ind*) <sup>2</sup> (*t*):

$$\begin{aligned} \mathbf{X}\_{1}^{(\text{ind})}(t) &= \int\_{0}^{t} dt' \mathbf{x}\_{1}^{(\text{ind})}(t') \\ &= \quad \mathbf{x}\_{1}(0) \sqrt{\frac{\pi t}{2c\_{1}}} e^{\frac{c\_{0}^{2}}{2c\_{1}}} \left[ \text{erf} \left( \frac{c\_{1}t - c\_{0}}{\sqrt{2c\_{1}}} \right) - \text{erf} \left( \frac{-c\_{0}}{\sqrt{2c\_{1}}} \right) \right], \\ \mathbf{X}\_{2}^{(\text{ind})}(t) &= \quad \int^{t} dt' \mathbf{x}\_{2}^{(\text{ind})}(t') \end{aligned} \tag{14}$$

$$\begin{aligned} \begin{pmatrix} \frac{\partial}{\partial t}^{\mathrm{(ind}}(t) &=& \int\_0^t dt' x\_2^{\mathrm{(ind}}(t')\\ &=& x\_2(0) \int\_0^t dt' e^{c'\_0 t' - c'\_1 t'^2/2 - c'\_2 t'^3/3} \end{pmatrix} \end{aligned} \tag{15}$$

.

where the error function is erf(*x*) = <sup>√</sup><sup>2</sup> *π x* <sup>0</sup> *dte*−*<sup>t</sup>* 2 .

To calculate *X*(*ind*) <sup>1</sup> (*t*), we used: *dt eat*−*bt*<sup>2</sup> = <sup>√</sup>*π<sup>e</sup> a*2 4*b* 2 <sup>√</sup>*<sup>b</sup>* erf-2*bt*−*a* 2 √*b* 

The model (11) was implemented using Python. We performed simulations taking the same scenarios adopted by Huang et al. and described in Section 3. Figure 1 shows the short-term population dynamics graphs of the predator–prey system—the number of organisms (No./m3) versus time (months)—where we observe two patterns: 1—the growth of the prey population at a rate *β*1(*t*) accompanied by the extinction of predators and 2—the extinction of both species.

As an example of the first pattern observed when *α* = *α* = 0, we consider the case where effects of MP particles on prey and predator are weak (*r*<sup>11</sup> = *r*<sup>21</sup> = 0.1). Since prey and predator are independent species and, according to the result of Huang et al., predators are more vulnerable than prey to the effects of MP particles, we have the growth of prey at a rate *β*1(*t*) accompanied by the expected extinction of predators (Figure 1a). This behavior was observed in two other scenarios: without the influence of MP particles (*C*<sup>1</sup> = *C*<sup>2</sup> = *r*<sup>11</sup> = *r*<sup>21</sup> = 0) and when predators have much larger response strength than prey to MP particles (*r*<sup>11</sup> = 0.01 and *r*<sup>21</sup> = 0.1), in which the effects of microplastics are not strong enough to stop the growth of the prey population. When the response intensities of MP particles on prey and predator are both increased to *r*<sup>11</sup> = *r*<sup>21</sup> = 5.0, we have an example of the second pattern, in which the effects of microplastics lead both species to extinction (Figure 1b). Prey population growth is substantially affected by this increase of *r*11, and the number of prey rises at *t* = 7.5, peaks at *t* = 13.7, and decreases to zero. All other simulated scenarios depict this complete extinction. In these cases, the increase of the adverse effect of reduced predatory performance (*d*2) and the consequent lost amount of

prey eaten by the predator (*d*3) led to an increasingly early extinction of prey, characterized by curves with increasingly less pronounced peaks.

**Figure 1.** Short-term population dynamics of the predator–prey for *α* = *α* = 0. Predator and prey have the same response strength to MP particles, i.e., Δ = *r*11/*r*<sup>21</sup> = 1.0. Abscissa measures time in months, and ordinate measures the number of organisms (No./m3). The populations of prey (*x*1) and predator (*x*2) are represented by the blue dotted line and the orange continuous one, respectively.

#### *3.2. Reduction of Predatory Performance with No Prey Eaten by Predator: α* = 0

Setting *α* = 0 leads to the second case where the equations of the system (10) decouple. In this way, we have

$$\begin{cases}
\dot{\mathbf{x}}\_1(t) = \mathbf{x}\_1(t)\boldsymbol{\beta}\_1(t) \\
\dot{\mathbf{x}}\_2(t) = \mathbf{x}\_2(t)[\boldsymbol{\beta}\_2(t) - a'\mathbf{x}\_1(t)]
\end{cases}.\tag{16}$$

In the original system (1), since *a*1*x*1*x*<sup>2</sup> refers to the lost amount of prey eaten by predators and *d*3*x*1*x*<sup>2</sup> refers to the decreasing of that amount, *α* = *a*<sup>1</sup> − *d*<sup>3</sup> = 0 describes the regime where the prey population does not decline due to predator feeding (a "total decreasing").

For *α* = 0, the population dynamics of prey are not affected by the predator population size. Nevertheless, the dynamics of predators are affected by the presence of prey. That happens because the model assumes a specific effect that affects the lost amount of prey eaten by predator (denoted by *d*3) and a distinct effect that affects the number of predators due to their reduced predatory performance (denoted by *d*2). Thus, even if there is not a lost amount of prey eaten by predators (*α* = *a*<sup>1</sup> − *d*<sup>3</sup> = 0), the adverse effect of reduced predatory performance may not be strong enough that the number of predators does not depend on the number of prey (*α* = *d*<sup>2</sup> − *a*2). In this scenario, the effects of microplastics make only the prey population independent.

In this scenario, since *x*1(*t*) = *x* (*ind*) <sup>1</sup> (*t*) = *<sup>x</sup>*1(0)exp *c*0*t* − *c*1*t* 2/2 , we have

$$\begin{aligned} \mathbf{x}\_2(t) &= \mathbf{x}\_2(0) \mathbf{e}^{\int\_0^t dt' \left[\beta\_2(t') - \mathbf{a}' \mathbf{x}\_1^{(ind)}(t')\right]} \\ &= \mathbf{x}\_2^{(ind)}(t) e^{-\mathbf{a}' \mathbf{X}\_1^{(ind)}(t)} \end{aligned}$$

where *x* (*ind*) <sup>2</sup> (*t*) = *<sup>x</sup>*2(0)exp *c* <sup>0</sup>*t* − *c* 1*t* 2/2 <sup>−</sup> *<sup>c</sup>* 2*t* 3/3 .

The model (16) was implemented using Python. We performed simulations taking the same scenarios adopted by Huang et al. and described in Section 3. Figure 2 shows the short-term population dynamics graphs of the predator–prey—the number of organisms (No./m3) versus time (months)—where we observe the same pattern, characterized by the growth of the prey population at a rate *β*1(*t*) accompanied by the predator "population explosion".

**Figure 2.** Short-term population dynamics of the predator–prey for *α* = 0. Predator and prey have the same response strength to MP particles, *r*<sup>11</sup> = *r*<sup>21</sup> = 0.1, i.e., Δ = *r*11/*r*<sup>21</sup> = 1.0. Abscissa measures time in months and ordinate measures the number of organisms (No./m3). The populations of prey (*x*1) and predator (*x*2) are represented by the blue dotted line and the orange continuous one, respectively.

On one hand, the prey population behaves in these scenarios in the same way as in each of the situations shown in Figure 1, since, in both cases, this population is independent of the predator population and the harmful effects of MP particles are the only factor limiting its growth. On the other hand, the predator population is substantially benefited by its dependence on prey. Figure 2 shows that, in the simulation range allowed by the Python compiler, the number of predators rises very quickly at around *t* = 0.4. This "explosive growth" was also observed in all other cases where *α* = 0, in which it was noted that it happens at around *t* = 0.5 when the values of *d*1, *d*2, and *d*<sup>3</sup> are increased.

*3.3. Reduction of Predatory Performance with No Increase in the Number of Predators Due to the Feeding of Prey: α* = 0

A third way to decouple the equations of the system (10) is to make *α* = 0:

$$\begin{cases}
\dot{\mathbf{x}}\_1(t) = \mathbf{x}\_1(t)[\beta\_1(t) - \alpha \mathbf{x}\_2(t)] \\
\dot{\mathbf{x}}\_2(t) = \mathbf{x}\_2(t)\beta\_2(t)
\end{cases}.\tag{17}$$

Since *a*2*x*1*x*<sup>2</sup> is the increasing number of predators due to the feeding of prey, and *d*2*x*1*x*<sup>2</sup> is the decreasing number of predators due to the adverse effect of reduced predatory performance, *α* = *d*<sup>2</sup> − *a*<sup>2</sup> = 0 defines a regime of reduction of predatory performance, where there is no increase in the number of predators due to the feeding of prey. Here, the exposure to MP particles makes only the predators population independent.

Since *x*2(*t*) = *x* (*ind*) <sup>2</sup> (*t*) = *<sup>x</sup>*2(0)exp *c* <sup>0</sup>*t* − *c* 1*t* 2/2 <sup>−</sup> *<sup>c</sup>* 2*t* 3/3 , we have

$$\begin{array}{rcl} \mathbf{x}\_1(t) &=& \mathbf{x}\_1(0)e^{\int\_0^t dt'\beta\_1(t') - \mathbf{a}\int\_0^t dt'\mathbf{x}\_2^{(ind)}(t')}\\ &=& \mathbf{x}\_1^{(ind)}(t)e^{-\mathbf{a}X\_2^{(ind)}(t)}. \end{array} \tag{18}$$

The model (17) was implemented using Python. We performed simulations taking the same scenarios adopted by Huang et al. and described in Section 3. Figure 3 shows the short-term population dynamics graphs of the predator–prey-number of organisms (No./m3) *versus* time (months)—where we observe the same two patterns of the Section 3.1, where *α* = *α* = 0.

**Figure 3.** Short-term population dynamics of the predator–prey for *α* = 0. Predator and prey have the same response strength to MP particles, i.e., Δ = *r*11/*r*<sup>21</sup> = 1.0. Abscissa measures time in months and ordinate measures the number of organisms (No./m3). The populations of prey (*x*1) and predator (*x*2) are represented by the blue dotted line and the orange continuous one, respectively.

When *α* = 0, the prey population is the only one affected by the presence of the other species, and predator extinction occurs early in every scenario, since there is no increase in its number due to the feeding of the prey, besides being under the harmful effects of the MP particles (Figure 3). Thus, the population dynamics of the system are similar to those in Figure 1. However, it can be seen in each of these situations that the dependence on the predator, even weak, led to an early extinction of the prey, which also presented curves with less pronounced peaks.

A possible next step for a system of ODEs such as (10) is to combine the two first-order differential equations into a single second-order one. To do so, we first write *x*2(*t*) as a function of *β*1(*t*) and the ratio *x*˙1/*x*<sup>1</sup> and differentiate *x*˙1(*t*) in Equation (10), leading to: *<sup>x</sup>*¨1(*t*) = *<sup>x</sup>*˙1(*t*)[*β*1(*t*) <sup>−</sup> *<sup>α</sup>x*2(*t*)] + *<sup>x</sup>*1(*t*)[*β*˙ <sup>1</sup>(*t*) − *αx*˙2(*t*)], and then to: *x*¨1(*t*) = *x*˙1(*t*)*β*1(*t*) + *x*1(*t*)*β*˙ <sup>1</sup>(*t*) − *α*{*x*˙1(*t*) + *x*1(*t*)[*β*2(*t*) − *α x*1(*t*)]}*x*2(*t*). Using the calculated *x*2(*t*), one obtains:

$$\ddot{\mathbf{x}}\_1(t) = \beta\_2 \dot{\mathbf{x}}\_1(t) + \left[\dot{\beta}\_1(t) - \beta\_1(t)\beta\_2(t)\right] \mathbf{x}\_1(t) + a'\beta\_1(t)\mathbf{x}\_1^2(t) - \left[-\frac{\dot{\mathbf{x}}\_1(t)}{\mathbf{x}(t)} + a'\mathbf{x}\_1(t)\right] \dot{\mathbf{x}}\_1(t), \quad t \in [0, T]$$

which does not depend on *α*. Analogously, if one had written the second-order equation for *x*2(*t*), it would be independent of *α* . This issue stresses that possibly Huang et al. should be reviewed.

#### **4. Discussion**

Rewriting the Huang et al. four equations model, reducing it to a two equations one and including time-varying intraspecies coefficients, allowed us to obtain an equivalent model with a format analogous to the standard LVM. This first step clarified the meanings of the model parameters and led us to explore and solve analytically special, and simpler, cases of ecological relationships. Bioaccumulation of microplastics through food web, and its consequent accelerated accumulation and magnification on predator [17–19], is a known effect that we analytically showed based on this model, calculating the time threshold from which their concentration and changing rate of the total amount are greater in predators than in prey.

The three situations that we studied are characterized by some possible behavioral abnormalities caused by the exposure to MP particles. More specifically, we considered effects of reduction of predatory performance [26], where prey is not affected by predator population size (*α* = 0), and/or there is no increase in the number of predators due to the feeding of prey (*α* = 0). In our study, theses effects are considered so strong that they can make one of the species independent, or even make both independent of each other. In the cases where *α* = *α* = 0 or only *α* = 0, results reveal basically two patterns, depending on the response strength to MP particles: 1—the growth of prey at a rate *β*1(*t*) accompanied by the extinction of predators, where the effects of microplastics are not strong enough to stop the growth of the prey population; and 2—the extinction of predator followed by the eventual extinction of prey. The most interesting, and counterintuitive, result of these three cases appears when *α* = 0, where there is predators "population explosion", showing the predator population is substantially benefited by its dependence on prey. We also point out that the Huang et al. model is not able to address the Allee effect [33,34]. It is important to emphasize that the underlying structure of the model is exponential, and the three special situations that we studied are limit cases, which lead to very different predictions from the original study and result in the extremely large numbers of organisms shown in Figures 1–3.

#### **5. Conclusions**

We rewrite the Huang et al. model by reducing its number of equations and defining new parameters, making it analogous to the standard LVM with time varying intraspecific coefficients. The time threshold from which the MP particles concentration and changing rate of its total amount are greater in predators than in prey was calculated. We solved analytically specific and simpler ecological situations where the effect of MP particles cause severe abnormal behavior on predator and prey, leading them to become independent of each other. Our simulations reveal a counterintuitive result when toxicological effects of MP particles cause a total interruption of prey feeding ability, which can produce a predator "population explosion".

Our advance in the analytical treatment of the Huang et al. model opens space for further refinement in the study of predator–prey models under toxicological effects of MP particles, either exploring the second-order equation that we propose or modifying the original model to further reduce its number of parameters, embedding in the time-varying intraspecies coefficients all the adverse effects caused by MP particles.

**Author Contributions:** A.S.M. obtained the analytical solutions, designed the research, wrote the manuscript and was responsible for supervision, project administration, and funding acquisition. M.M.R.S. introduced the Huang et al. paper to the research group, and wrote and validated the first version of the numerical code. L.M.D.C. introduced the Huang et al. paper to the research group, wrote, and validated the first version of the numerical code. J.R.A. verified analytical solutions, identified the critical values of Equations (7) and (8), and wrote the manuscript. L.S.D.S. verified analytical solutions; wrote a second and independent version of the numerical code, validated it and understood the numerical overflow relative to Figure 2, interpreted the model parameters with ecological regimes, and wrote the first version of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** José Renato Alcarás and Mateus Mendonça Ramos Simões have been funded by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brasil (CAPES)—Finance Code 001. José Renato Alcarás acknowledges support of Fundação de Amparo à Pesquisa do Estado de São Paulo-Brazil (FAPESP) (Grant No. 2018/21694-4). Mateus Mendonça Ramos Simões acknowledges support of Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brazil (CAPES) (Grant No. 88887.570033/2020-00) Alexandre Souto Martinez acknowledges Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (Grant No. 309851/2018-1) and NAP-FisMed for support.

**Data Availability Statement:** The numerical code in Python is available on demand.

**Acknowledgments:** The authors thank Stephany Silva Xesquevixos de Siqueira for fruitful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

MP Microplastic

LVM Lotka–Volterra model

#### **References**

