Circuit Design and Experimental Investigations for a Predator–Prey Model

Consider the predator–prey model with Allee Effect described by Ben Saad and Boubaker, (2015, 2017): (1) {dx1dt=x1(x1−ℓ)(k−x1)−x1x2,dx2dt=e(x1−m)x2, where x₁ is the size of the prey population and x₂ the size of the predator population, ℓ is the Allee Effect threshold, k is the carrying capacity, e is the feeding efficiency of Lotka–Volterra model and m is the predator mortality rate. Let consider also the biologically meaningful conditions x₂ ≥ 0 and x₂ ≥ 0.

Two case studies for the Allee effect ℓ can be considered (Dhooge et al., 2008):

The Strong Allee effect when ℓ ∈ [0 1],

In the following, only the Weak Allee effect case study will be considered and the parameters k, e, and m are taken as k = e = 1 and m ϵ [0 1].

For k = e = 1, ℓ ∈ [−1 0] and m ϵ [0 1], the system (1) admits three equilibriums E_i(i = 1,…,3):

The zero equilibrium E₀ = (0, 0),

For i = 1, …, 3, the corresponding Jacobian matrix J_i and eigenvalues (λ_1i, λ_2i) are given by (2) J0=(−ℓ00−m)and{λ10=−l,λ20=−m (3) J1=((ℓ−1)−10(1−m))and{λ11=l−1,λ21=1−m (4) J3=((−2m2+m(ℓ+1))−m(m−ℓ)(1−m)0)and{λ12=tr(J3)+iΔ2λ22=tr(J3)−iΔ2 with tr(J₂) = m(ℓ+1 − 2m); Δ = [m(ℓ+1 − 2m)]²−4[m(m − ℓ)(1 − m)].

Based on the last results, the singularities and the phase portraits of the neighborhood of all equilibriums are summarized in Table 1.

The global dynamics of the system (1), already detailed in Table 1, will be proved numerically in this section using the Matcont software (Dhooge et al., 2008). Since the two first equilibriums E₀ and E₁ exist ∀m ∈ [0 1], only the singularity of the equilibrium E₃ will be shown numerically in the following for the initial condition (x₁ = 0.4; x₂ = 0.3). For m = −ℓ =0.2, Figure 1 shows the temporal evolution of the dynamics of system (1) whereas the phase portrait of the system is shown in Figure 2.

As it is shown, a heteroclinic orbit appears indicating the existence of a global bifurcation described by an attractor which bifurcates from an unstable focus equilibrium point to a heteroclinic cycle. The two populations’ dynamics oscillate between 0 and 1 density values. Figures 3 and 4 show the dynamic behavior of the system at the non-isolated point E₃ (0.4, 0.36) when m = 0.4.

For the second case study, the obtained behavior of the two populations is periodic which is described as a circle around the equilibrium E₃ (0.4, 0.36). When m = 0.6, the dynamic behavior of the system is shown in Figures 5 and 6.

The obtained phase portrait shows a stable focus. This latter closes to a fixed point with the following coordinates (0.6, 0.32) proven with the temporal evolution.

Fig. 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Fig. 6

The aim of this section is to design an analog circuit that can build the nonlinear terms according to system (1). Therefore, equivalent electronic circuits of these nonlinear terms are designed and simulated via MultiSIM software (Yujun et al., 2010). Due to the complexity of the equivalent circuit design, we designed first the equivalent circuit of each nonlinear term x² and x³. The first nonlinear term circuit needs one Multiplier AD6333. The multiplying core of this multiplier comprises a buried Zener reference providing an overall scale factor of 10 V. For that, an amplification of the output with a gain of 10 is required. Thus, the equivalent electrical model of x² is realized as follows: S1=x2=R210R1x2.

As shown in Figure 7, by considering the components datasheet, the electronic circuit of the nonlinear term x² is modeled by one Analog multiplier AD633AN, one Amplifier AOP, two capacitors (C₁, C₂) and two resistors (R₂, R₂). Due to the existence of internal loss in the electrical components, the ideal value of R₂ chosen to obtain the best result is 0.9 kΩ. However, for the nonlinear term x³, the electronic circuit is modeled in Figure 8 by two multipliers AD633AN, 2 Amplifiers AOP, four capacitors and four resistors, since its electrical model is realized as follows: S2=x3=R410R3S∗x=R4R2100R1R3x3.

The input which will be used for the equivalent circuit of system (1) is continuous. However, in this part, in order to verify the efficiency of the two nonlinear terms equivalent circuits, an alternative signal is used as an input with a weak frequency equal to 1 Hz and amplitude of 2 V. The square and the cube of the signal are presented with a pink curve in Figure 9 and Figure 10, respectively. For high frequencies (1 kHz, 10 kHz, 100 kHz), we obtained the same following results.

As it is shown, we have chosen the same scale for the two nonlinearities. For the first nonlinearity x², the obtained signal amplitude is equal to 4,108 V which obviously the square of the input signal amplitude. In addition, for the second nonlinearity x³, the result proves that the obtained signal is the cubes of the input signal since the obtained amplitude is equal to 7,980 V which is almost equal to 8 V the cube of the input amplitude.

Figures 11 and 12 show an experimental implementation for the last simulation results. The oscilloscope traces of the proposed circuit are illustrated in Figures 13 and 14.

Comparing the different results, it is proven that there is a good qualitative agreement between the numerical simulations with Matlab, the electrical simulations with MultiSIM and the experimental results for the two nonlinearities.

Fig. 7

Fig. 8

Fig. 9

Fig. 10

Fig. 11

Fig. 12

Fig. 13

Fig. 14

In this section, the agreement between biological theory and electronic experiments of the predator–prey model will be analyzed by considering the three cases studies presented numerically by MATCONT. Therefore, a transformation of the biological predator–prey model (1) to an equivalent electrical model is realized as follows: (5) {dx1dt=−1C1R1x13+1C1R2x12−1C1R3x1−1C1R4x1x2,dx2dt=1C2R5x1x2−1C2R6x2, where 1C1R1=1C1R4=1C2R5=1 ; 1C1R2=(1+l) ; 1C1R3=l ; 1C2R6=m . C₁ and C₂ are used for the integration of the circuit outputs in order to obtain as output the populations’ density x₁ and x₂.

The electronic circuits corresponding to these cases are designed by the software MultiSIM and presented in the following. We have used three multipliers, five AOP, two capacitors and 12 resistors. The resistors (R₁, R₂, …, R₆) and capacitors values are fixed with respect to the parameters values. The value of the two capacitors C₁ and C₂ is fixed at 100 nF. In addition, two interrupters S1 and S2 are used to introduce the initial conditions of the prey and predator density. In order to analyze the three case studies, we vary the value of the resistor R₆ which corresponds to the parameter m since R6=R5m . In the first case study, we consider the mortality rate of the predator m = 0.2. Therefore, the fixed value of R₆ in this case is equal to 50 kΩ. Figure 15 describes the electrical circuit of the predator–prey model (5).

Then, the temporal evolution and the phase portrait when m = 0.2 are presented in Figure 16 and Figure 17, respectively.

We obtained two phase-shifted alternative signals. The pink curve corresponds to the dynamic evolution of the prey while the blue one corresponds to the dynamic evolution of the predator. Based on the scale chosen in the temporal evolution, the amplitude of the blue and pink curves are almost equal to 0.956 V and 0.836 V, respectively. Comparing with the numerical simulation presented in Figures 1 and 2, we conclude that the permanent regime of the system dynamic obtained via MultiSIM is similar to that obtained using the Matcont software.

In the second case study, we consider the mortality rate of the predator m = 0.4. Thus, Figure 18 describes the electrical circuit of the model with R₆ = 25 kΩ. Then, the temporal evolution and the phase portrait are illustrated in Figure 19 and Figure 20, respectively.

As it is illustrated in the temporal evolution, the maximum amplitude of the prey population is equal to 0.488 V and that of the predator population is equal to 0.441 V. The obtained values are obviously very close to the numerical values shown in Figure 3. Furthermore, the phase portrait proves the existence of the center singularity which is demonstrated numerically via the Matcont software in Figure 4.

In the third case study, the mortality rate of the predator m is fixed at 0.6. Therefore, Figure 21 presents the electrical circuit with R₆ = 17 kΩ. Figures 22 and 23 illustrate the corresponding temporal evolution and phase portrait.

For the third case study, we change slightly the initial conditions (x₁ = 0.4, x₂ = 0.1) to obtain the clearest results. The obtained temporal evolution and phase portrait prove that the model dynamic tends toward a fixed point. In Figure 22, it is noted that the dynamic behavior oscillates during the transitional regime, whereas in the permanent regime, it becomes continuous reaching a constant value equal to 0.586 V for the prey population and to 0.326 V for the predator one. The pink curve corresponding to the prey population dynamic reaches a pick of 0.978 V. However, the blue curve corresponding to the predator population reaches a pick equal to x² = 0.399 V. This temporal behavior is described by a stable focus in the obtained phase portrait. These electrical results are almost similar to the numerical ones presented in Figures 5 and 6.

For the three cases study, we conclude that the electrical results obtained by the software MultiSIM obviously prove the numerical results obtained by the Matcont software.

Fig. 15

Fig. 16

Fig. 17

Fig. 18

Fig. 19

Fig. 20

Fig. 21

Fig. 22

In this section, an experimental implementation of the equivalent electrical circuit will be realized by placing the different electrical components on a bread board within the laboratory.

As it is mentioned previously, mathematical predator–prey model has several nonlinearities which are modeled by multipliers AD633 in the equivalent electrical circuit. The cascade of components AD633 and AOP with gain of 10 amplifies most likely the imperfection and the uncertainties of the electronic components and consequently distorts the experimental results. Therefore, in order to resolve this problem, the experimental validation is realized by using the technology STM3278 shown in Figure 24.

The experimental results of the three cases studies are obtained by using the oscilloscope and presented in Figures 25–30.

Comparing experimental results with numerical results obtained via the MATCONT software shown in Figures 1–6, we can conclude that good agreement is obtained between simulation results and experiments.

Fig. 23

Fig. 24

Fig. 25

Fig. 26

Fig. 27

Fig. 28

Fig. 29

Fig. 30

In this paper, an electronic circuit is designed for the model of a prey–predator model and its complex behavior is proved via numerical and electrical results. Furthermore, some experimental investigations are also shown and demonstrate that they can be used to characterize the ecological dynamics faster. In the future, we will try to prove, via experimental results, the chaotic dynamics of such nonlinear systems in the presence of seasonally effects in order to use such circuits for encryption/decryption fields.