Modulation of the mechanical properties of ventricular extracellular matrix hydrogels with a carbodiimide crosslinker and investigation of their cellular compatibility

1.   Introduction

In agronomy, tritrophic interactions between crop, herbivores and their natural enemies are one of the drivers of the crop yield. Understanding and manipulating these interactions in order to produce food more sustainably is the basic principle of biological control of pest. Plants continuously emit a blend of different Volatile Organic Compounds (VOCs). Among them, some are induced or produced in a higher quantity when the plant is attacked by an herbivore ^[1,2,3,4]. They may be considered as an indirect defence mechanism of the plant as they attract natural enemies toward the attacked plants and therefore to the pests themselves ^[1,2,5,6]. This ecological process is widely common since the range of natural enemies reacting to different plants VOCs includes many diverse predators, such as insect and mites, and many diverse parasites such as parasitoid wasps or flies, and even parasitoid nematods ^[1]. One specific example of pest-natural enemy is the interaction between aphids and ladybirds. For example, Coccinella septempunctata has been shown to be attracted by VOCs emitted by (ⅰ) barley Hordeum vulgare when infected by the bird cherry oat aphid Rhopalosiphum padi ^[7] and (ⅱ) several cruciferous vegetables Brassica spp. when infected by the peach potato aphid Myzus persicae ^[8].

Mathematical modelling of tritrophic interactions may be dated back to the seventies when three–species food chains were studied by means of ordinary differential equations (ODEs) as extensions of two–species interactions of predator–prey type ^[9,10]. Later on, in a very famous paper that inspired a wealth of subsequent research, Hastings and Powell ^[13] used such kind of models to describe the interaction between a resource, a consumer and a predator. Their model (say, HP model) was an evolution of the very well known two–dimensional Rosenzweig–MacArthur model ^[11]. An important outcome of the HP model is that nonlinear responses (in particular a Holling type Ⅱ (HTII) functional response was used) in the long–term behavior may favor the onset of chaotic dynamics. This result was somehow criticized by McCann and Yodzis ^[12] who argued that the parameter values used in ^[13] were `rather improbable' biologically. However, McCann and Yodzis also indicated new parameters that also produce chaos. Therefore, a general conclusion gained from these studies was that even relatively simple ODE food chains model may show a very rich dynamics. These models are able, in particular, to describe natural food webs chaotic dynamics that `ought to be quite common when the resource productivity is sufficiently high' ^[12]. Since then, chaos in multi–species systems has been the subject of many investigations (see e.g. the extensive review of chaotic models in ecological systems contained in ^[14]). The HP model itself has been studied through bifurcation theory ^[15,16] and numerical simulations ^[17]. On the other hand, the HP model has been later revisited to represent many different scenarios. This was done, for example, by Maiti and co-workers ^[18] who studied the interaction of tea plant, pest and predator. Pal and co-workers ^[19] used that model to prove that omnivory in tritrophic food chain stabilizes the system from chaos to order. Zhang and Georgescu ^[20] used an HP–like model coupling the within-pest virus dynamics with the crop-pest virus dynamics. Gakkar and Singh ^[21] extended the HP model to a four dimensional model by introducing an additional predator in the top prey. Through numerical simulations they showed that for a given set of parameter values the modified food web model exhibits stable dynamics in contrast to chaotic dynamics occurring in the three species HP model. The HP model with more general functional responses has also been studied in both deterministic ^[22,15] and stochastic enviroments ^[23].

A special application of the tritrophic models is concerned with the effects of chemicals released by populations involved in the tritrophic interaction. Mandal and co-workers considered the release of toxins in an aquatic ecosystem ^[24], while Liu and co-workers proposed a model where allelochemicals are released by plant species as a protection against herbivores or pests ^[25]. Later on, this model has been revisited by Fergola and Wang ^[26].

In this paper, we consider a model given by three ordinary differential equations describing the tritrophic interaction between crop, pest and the pest natural enemy, in which the release of Volatile Organic Compounds (VOCs) by crop to attract the pest natural enemy is explicitly taken into account. The model we propose is an extension of the HP model ^[13], in the sense that an extra nonlinear term is added to the equation ruling the dynamics of the highest trophic level to describe the effects of VOCs. On the other hand, it is also a more complex variant of the plant–herbivores–enemies considered in ^[26,25], where the interspecific interactions are described by bilinear mass–action terms.

The rest of the paper is organized as follows: the model is introduced in Section 2. In Section 3 we present the equilibria analysis for feasibility and local stability. Section 4 is devoted to numerical simulations to assess the VOCs effects on the model dynamics. Conclusions are given in Section 5.

2.   The model

We assume that the size of the interacting populations, i.e. the crop, the pest/aphid and the pest natural enemy, are represented by the time-dependent functions $x(t)$, $y(t)$ and $z(t)$, respectively. As possible example of a concrete tritrophic interaction including VOCs we consider barley, H. vulgare ($x$), the aphid R. padi ($y$) and the ladybird C. septempunctata ($z$).

Unlike the model proposed in ^[26] we assume that the production of volatiles by the plants in the unit time is limited. To this end, we replace the mass action terms used in ^[26] by HTII responses. As a consequence, the dynamics is ruled by the following three ODEs, where all the parameters, whose meaning and baseline values are summarized in Table 1, are positive except $b \geq 0$ and $c \geq 0$:

Table 1.  Model variables and parameters: notation, biological description and baseline values used in the simulations.

Symbol	Description	Baseline value
$x$	crop population size	$0.2$ initial value
$y$	aphid population size	$0.2$ initial value
$z$	aphid-natural enemy population size	$0.2$ initial value
$r$	crop growth rate	$0.1$
$K$	crop carrying capacity	$1$
$a$	maximal harvesting rate of crop by aphids	$0.1$
$e$	crop to aphids conversion (yield)	$0.4$
$m$	aphids' natural mortality rate	$0.01$
$p$	maximal uptake rate of aphid by aphid-natural enemy	$0.01$
$h$, $k$, $\ell$	half saturation constants	$0.5$
$b$	attraction constant due to VOCs	$0.25$
$c$	enhanced attraction rate of aphid-natural enemy by VOCs released by crops under aphid attack	0.5
$q$	aphids to aphid-natural enemy conversion (yield)	0.5
$n$	aphid-natural enemy mortality rate	0.3

---

 | Show Table

DownLoad: CSV

In (2.1) the first equation rules the crop dynamics. The crop is assumed to grow logistically with net reproduction rate $r$ and carrying capacity $K$. The crop is attacked by the aphids population, an action that is modeled in the last term of the first equation by means of an HTII response function, to account for feeding satiation of the insects. The dynamics of the aphid is ruled by the second equation, where the conversion coefficient $e$ is introduced in the uptake of the consumed crop. Further, the aphids die with natural mortality rate $m$; they are also attacked by their natural enemies, i.e. the C. septempunctata population, at rate $p$, again using an HTII function. Finally, the third equation describes the dynamics of pest natural enemy. The first two terms represent, respectively, the pest enemy attraction due to: i) VOCs basic plant release, at rate $b$; ii) VOCs plant release due to pest attack, at rate $c$. As mentioned above, the latter process is modelled by an HTII function. Pest-enemy reproduces by feeding on the aphids, with conversion factor $q$, and dies with mortality rate $n$. Note that the input of new pest natural enemies, expressed by the first two terms in the last equation, implies the existence of an unlimited external reservoir of aphid-natural-enemies which is not explicitly modelled in (2.1).

A conceptual diagram of the model structure is shown in Figure 1, where the state variables and the related processes are represented.

Figure 1.  Schematic representation of the inter-specific interactions between plant, herbivore pest and pest-natural enemy. The solid lines represent the feeding processes (1, 2), the dotted lines represent the increase of volatile compounds due to: basic plant release (3) and plant release during pest attack (4). The volatile compounds attract the pest-natural enemy from the external reservoir, which is not explicitly modeled in the dynamical system (2.1), whose state variables $x$, $y$ and $z$ are represented by the rectangular boxes. Note, however, that volatiles are not explicitly modeled as system variables.

DownLoad: Full-Size Img PowerPoint

3.   Equilibria and stability

Possible equilibria are only the ecosystem collapse $E_0 = (0, 0, 0)$, the aphid-free point $E_1 = (x_1, 0, z_1)$ and coexistence $E_* = (x_*, y_*, z_*)$. In particular

It follows that the point $E_1$ always exists. For coexistence we can find some sufficient conditions for its existence as follows.

We find $y$ from the first equilibrium equation,

Note that this is a quadratic function in $x$, which is nonnegative only for $0\le x \le K$. Solving directly for $z$, from the second equilibrium equation we find

This function is nonnegative for $aex \ge m(h+x)$. Let $\bar x$ be the unique solution of the equation

The conditions for which $z\ge 0$ are therefore:

We need to take an interval in which both $y\ge 0$ and $z\ge 0$. Therefore, for (3.3) we need $\bar x \le K$, while for (3.4) it follows $\bar x \le 0$ and no feasible interval can arise. Imposing that $\bar x$ lies in the interval of interest $[0, K]$ we obtain from (3.3)

Define now a new function $f(x)$ by the right hand side of the third equilibrium equation in which both $y$ and $z$ are substituted with the above expressions (3.1) and (3.2). Explicitly,

It is then easily checked that

in view of the fact that $\bar x \le K$. On the other hand,

so that $f(K) \le 0$ if we impose

and by its continuity, the function $f$ must have a zero $x_*$ in the interval $[0, K]$. The remaining coordinates of $E_*$ are obtained then from (3.1) and (3.2). Thus sufficient conditions for the existence of the equilibrium at which all populations coexist are (3.5) and (3.7).

For stability we need the Jacobian of (2.1):

Evaluation at $E_0$ shows that this point is unconditionally unstable, in view of the positive eigenvalue $r > 0$. At $E_1$ instead, two eigenvalues are negative, $-r$ and $-n$, while from the third one, the condition for local stability is obtained,

Comparing (3.7) and (3.9), if $b = 1$ there would be a transcritical bifurcation for which a coexistence equilibrium would arise when $E_1$ would become unstable.

The local stability condition of $E_1$, (3.9), shows that the disappearance of aphids is enhanced by their larger natural mortality $m$, by a more effective basic release of volatiles $b$, lower C. septempunctata's natural mortality $n$ and aphids attack rate on crops $a$.

The role of the crops size, expressed by their carrying capacity is not clear and therefore we turn to simulations to try to elucidate it.

In Figure 2 we show the equilibria and oscillations in the $e-K$ parameter space. It can be seen that the aphid population R. padi vanishes in a large portion of the parameter space, except for the lower values of the conversion coefficient $e$ and with growing values of the crops carrying capacity $K$. In this range coexistence of the whole ecosystem occurs through sustained oscillations. Their amplitudes appear to increase with increasing values of the conversion coefficient $e$. Here we used the set of initial conditions $x(0) = 1.7$, $y(0) = 2.8$, $z(0) = 3.0$, but another run, not reported, with $x(0) = 0.2$, $y(0) = 0.2$, $z(0) = 0.2$ (which we will use in the next section) does not lead to any significant change in the plots of these surfaces. Also, we see that the C. septempunctata population tends to grow with the size of the crops $K$. This of course assumes that this species is not a specialist predator on R. padi, because it thrives also when the pest disappears. In nature, in Europe C. septempunctata can live on other aphid species that also feed on cereals such as Sitobion avenae and Methopolophium dirhodum, ^[27]. To achieve eradication the upper portion of the parameter space must therefore be used.

Figure 2.  Equilibria and oscillations in the $e-K$ parameter space. The other parameter values are given in Table 1. Initial conditions: $x(0) = 1.7$, $y(0) = 2.8$, $z(0) = 3.0$. The single surface denotes the value of the population at equilibrium. When two surfaces appear, they represent the minimum and maximum values of the population persistent oscillations.

DownLoad: Full-Size Img PowerPoint

4.   Numerical Simulations

We perform a few further simulations to study the model behavior. Unfortunately, there is a lack of full data sets representing our specific biological system. Therefore in order to highlight the various possible dynamics shown by model (2.1) in the numerical simulations we use hypothetical parameter values, as has been previously done in similar investigations ^[26,13,9]. The biological meaning of each parameter can be obtained from the model, as reported in Table 1. Even though the model is theoretical and not yet related to field data, a concrete example of this ecosystem is represented, as mentioned in Section 2, by the H. vulgare, R. padi and C. septempunctata interactions. In this case, since the H. vulgare is an annual plant, the system (2.1) models only the vegetative growth phase of the plant. Indeed, recall that in the absence of interaction with the pest, the crop growth is logistic.

The first investigation deals with the effects of the conversion coefficient $e$ appearing in the aphids equation (the second equation of (2.1)) on the model dynamics. It can be seen that increasing the value of $e$ destabilizes the model and oscillations arise as shown in Figures 3 and 4. The threshold value in this case is approximately $e_c = 0.42$. A further increase of $e$ causes an increase in the amplitude of the oscillations, Figure 5.

Figure 3.  Dynamics of model (2.1) for $e = 0.4$. The other parameter values and initial conditions are given in Table 1. The state variables approach a steady state.

DownLoad: Full-Size Img PowerPoint

Figure 4.  Dynamics of model (2.1) for $e = 0.42$. The other parameter values and initial conditions are given in Table 1. Oscillations occur.

DownLoad: Full-Size Img PowerPoint

Figure 5.  Dynamics of model (2.1) for $e = 0.5$. The other parameter values and initial conditions are given in Table 1. Compared with the dynamics shown in Figure 4, here oscillations increase their amplitudes.

DownLoad: Full-Size Img PowerPoint

When $b = c = 0$, i.e. no pest natural enemies are attracted either by basic crops release or by VOCs enhancement due to the presence of the aphids attacking the crops. The natural-pest enemy rapidly disappears and an oscillatory dynamics involving just crops and aphids arises, Figure 6. Note that in this case the model reduces to the classical HP model.

Figure 6.  Dynamics of model (2.1) for $b = c = 0$. The other parameter values and initial conditions are given in Table 1. The model becomes the classical HP model in the regime of sustained oscillations. The natural-pest enemy rapidly disappear and an oscillatory dynamics involving only barley crops and aphids arises.

DownLoad: Full-Size Img PowerPoint

Then, we explore the effects due to the increase of the enhanced volatile release rate $c$ when crops are under aphids attack in the ladybird equation (the third equation of (2.1)). This is shown in Figures 6–8. When $c = 0.25$ the system shows an oscillatory behavior involving all the three state variables (Figure 7), while a coexistence equilibrium of all the three interacting populations can also be achieved as it can be observed for $c = 0.44$ (Figure 8).

Figure 7.  Dynamics of model (2.1) for c = 0.25. The other parameter values and initial conditions are given in Table 1. Oscillatory behavior involving all the state variables.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Dynamics of model (2.1) for $c = 0.44$. The other parameter values and initial conditions are given in Table 1. The state variables approach a coexistence equilibrium.

DownLoad: Full-Size Img PowerPoint

Finally, in order to analyze the separate effects of basic plant release and plant release during the pest attack, we perform different simulations by varying the values of the parameters $b$ and $c$. The simulations show that in the absence of release of volatiles by the crops (i.e. $b = 0$) the solution is a periodic solution of period one, see Figure 9 left panel. Similarly, when $c = 0$, the solution is qualitatively similar to the previous case with (possibly stable) oscillations, see Figure 9 right panel. Instead, considering both effects ($b\neq0$ and $c\neq0$) the populations achieve a coexistence equilibrium, see Figures 8 and 10.

Figure 9.  The figure shows the phase space of (2.1). Left panel: $b = 0$; right panel: $c = 0$. The other parameter values and initial conditions are given in Table 1. In both cases, period one oscillations are observed.

DownLoad: Full-Size Img PowerPoint

Figure 10.  Phase space of (2.1) for b = 0.25, and c = 0.44, while the other parameter values and initial conditions are given in Table 1.

DownLoad: Full-Size Img PowerPoint

5.   Conclusions

In this paper we have introduced a mathematical model to describe the tritrophic interaction between crop, pest and the pest natural enemy, in which the release of Volatile Organic Compounds (VOCs) by crop is explicitly taken into account. The VOCs have the role of attracting the pest natural enemy. Our model fits in the model scheme proposed by Hastings and Powell ^[13] and is a particular variant of the model proposed by Fergola and Wang ^[26].

Specifically, compared to ^[26], there are here two novel features. First, we consider saturation in feeding by the two insects, population pest and pest-natural enemy. Second, although recruitment of new pest-natural enemies is already present in ^[26], we assume that it is enhanced by the presence of aphids on crop. As a consequence, the modification introduced in the first term of the last equation entails that the crop-only equilibrium found in ^[26] is not possible here anymore.

The dynamics of the model is then clear, as only either the aphids-free equilibrium $E_1$ is allowed, when (3.9) holds, or coexistence of the whole ecosystem. In addition, the latter is seen to bifurcate, for instance by an increase of the conversion factor $e$ of crop into new aphids. These oscillations may drive the population close to zero at certain times, see e.g. the central frame of Figure 2 for the smaller values of $e$. In this way thus, theoretically, the oscillations may in principle threaten the survival of the ecosystem. Indeed possible stochastic environmental perturbations may change the model parameter values end entail the vanishing of some ecosystem population. From the biological point of view this is unlikely, though, because the colonization of the crop by the aphids is a phenomenon that occurs through the growing season of the crop.

The analysis of the local stability condition of the $y$–free equilibrium $E_1$ shows that several factors may enhance the disappearance of aphids as the larger natural mortality: a more effective release of volatiles and lower C. septempunctata's natural mortality and/or aphids attack rates on crops. Also, we have shown that the C. septempunctata population tends to grow with the crops carrying capacity $K$. Our analysis indicates that increasing the conversion coefficient $e$ may destabilize the model dynamics with the onset of sustained oscillations and that to ensure pest eradication a relatively high value of $e$ should be sought. However this is a rather difficult situation to achieve, as this is an intrinsic parameter of the aphid population. It is therefore better to act on the remaining model parameters appearing in the condition (3.9) to try to combat the pest.

Note further that C. septempunctata thrives also in the absence of the aphids R. padi because although it may lack the food explicitly modeled in the system, there is a continuous replacement by new incoming entries, attracted by the volatiles released by the plants, a fact that is expressed in the model by the first recruitment term in the last equation.

As far as the role of volatiles by crop is concerned, VOCs have been shown to possess a beneficial effect on the environment since their release may be able to stabilize the model dynamics. In particular, this effect may occur, in the parameter range that has been explored, only when both the phenomena of basic volatile release and volatile release during the pest attack are present, i.e. when $b\neq0$ and $c\neq0$.

Although the analysis and numerical simulations provide a satisfactory description of model dynamics, we are aware that this is a preliminary study. First of all, we have used different sets of hypothetical data. However, only field data, when available, will clarify the `right' dynamics that can be predicted from this model. From a mathematical point of view, the nonlinearities make difficult to obtain a complete qualitative analysis. The nature of the oscillations, the occurrence of Hopf bifurcation and the parameter values range ensuring possible chaotic dynamics will be the subject of further investigations.

Acknowledgements

This research is performed under the auspices of COST Action FA1405 "Using three-way interactions between plants, microbes and arthropods to enhance crop protection and production". FG and EV are members of the INdAM Research group GNCS, BB is member of the INdAM Research group GNFM. SS was supported by the Research Council of Norway through the SMARTCROP project (project number 244526). The authors warmly thank all the three referees for their comments that helped in improving the paper presentation.

Conflict of interest

The authors declare there is no conflict of interest.