Proof. The results follow from inspections of the eigenvalues of the matrix JE2 and Theorem 2 (see [16, 18]). ∎

Definition 4. [17] The flow F will be called uniformly persistent if there exists ε0 > 0 such that

The following theorem shows the existence of the equilibrium point E2 using persistence

sxy.

limt!þ∞ inf s tð Þ > k, limt!þ∞ inf x tð Þ > k, limt!þ∞ inf y tð Þ > k,

Proof.

Theorem 4. Suppose the assumptions such as

We can see from both Jacobian matrices JE1 and JE2 that there are no parameter b involved when the predator organisms y is zero. This means that intratrophic predation does not affect the local stability and existence of E1 and E2 when no predator organisms involved.  Theorem 5. If

$$\frac{f\_1'(\zeta\_\*)(\zeta\_\*-1) + f\_1(\zeta\_\*) - D\_1}{D\_1} < 0 \text{ and } \frac{f\_2\left(\frac{1-\zeta\_\*}{D\_1}\right)(1+2D\_1-\zeta\_\*-D\_1\zeta\_\*)-1}{D\_1-\zeta\_\*+1} < 0,$$

then the equilibrium point E2 ζs; ð Þ 1�ζs D1 ; 0 is a hyperbolic saddle and is repels locally in y-direction. Particularly, the dimension of the stable manifold Ws ð Þ E2 and the unstable manifold Wuð Þ E2 are given by

$$
\dim \mathcal{W}^s(E\_2) \quad = 2 \text{ and } \dim \mathcal{W}^u(E\_2) = 1.
$$

Proof. The results follow from inspections of the eigenvalues of the matrix JE2 and Theorem 2 (see [16, 18]). ∎

Definition 4. [17] The flow F will be called uniformly persistent if there exists ε0 > 0 such that for all x∈ E0 , limt!∞ dð Þ πð Þ x; t ; ∂E ≥ ε0:

The following theorem shows the existence of the equilibrium point E2 using persistence analysis.

Theorem 6. Assume that

i. Lemma 1 being holds,

ii. E2 is a unique hyperbolic saddle point in Rþ sxy and repels locally in y-direction (as in Theorem 5),

iii. no existence of periodic orbits in the planes of Rþ sxy.

Then

$$\lim\_{t \to +\infty} \inf s(t) > k, \lim\_{t \to +\infty} \inf x(t) > k, \lim\_{t \to +\infty} \inf y(t) > k.$$

where k > 0:

Particularly, the system (3) exhibits uniform persistence and the equilibrium point E2 ζs; ð Þ 1�ζs D1 ; 0 exists.

Proof. The result follows from the Definition 4, which defines uniform persistence by Butler et al. [17] and Nani and Freedman [18]. ∎ (iii) Stability analysis of E3 s∗; x∗; y∗ ð Þ:

#### The Jacobian matrix at E3 is

$$J\_{E\_{3}} = \begin{bmatrix} -f\_{1}(s^{\*})x^{\*} - 1 & -f\_{1}(s^{\*}) \\ f\_{1}(s^{\*})x^{\*} & f\_{1}(s^{\*}) \\ D\_{1} & D\_{1} - f\_{2}(x^{\*})y^{\*} \left[ \left(\frac{D\_{1}}{1+x^{\*}+D\_{1}by^{\*}}\right) - \left(\frac{D\_{1}x^{\*}}{(1+x^{\*}+D\_{1}by^{\*})^{2}}\right) \right] - f\_{2}(x^{\*})y^{\*} \left(\frac{D\_{1}x^{\*}}{1+x^{\*}+D\_{1}by^{\*}} + 1\right) - 1 \\ 0 & f\_{2}(x^{\*})y^{\*} \left[ \left(\frac{D\_{1}}{1+x^{\*}+D\_{1}by^{\*}}\right) - \left(\frac{D\_{1}x^{\*}}{(1+x^{\*}+D\_{1}by^{\*})^{2}}\right) \right] + f\_{2}(x^{\*})y^{\*} \left(\frac{D\_{1}x^{\*}}{1+x^{\*}+D\_{1}by^{\*}} + 1\right) \\ 0 & 0 \\ \left(\frac{f\_{2}(x^{\*})D\_{1}^{2}bx^{\*}y^{\*}}{(1+x^{\*}+D\_{1}by^{\*})^{2}}\right) - f\_{2}(x^{\*}) \left(\frac{D\_{1}x^{\*}}{1+x^{\*}+D\_{1}by^{\*}} + 1\right) \\ f\_{2}(x^{\*}) \left(\frac{D\_{1}x^{\*}}{1+x^{\*}+D\_{1}by^{\*}} + 1\right) - \left(\frac{f\_{2}(x^{\*})D\_{1}^{2}bx^{\*}y^{\*}}{(1+x^{\*}+D\_{1}by^{\*})^{2}}\right) - 1 \end{bmatrix}.$$

The eigenvalues of JE3 are resulted from the characteristic equation below

$$c\lambda^3 + c\_4\lambda^2 + c\_5\lambda + c\_6 = 0,\tag{6}$$

λ8 ¼ � 1 2D1

0

BBBBBB@

þ

þf 0 1 s ∗ ð ÞD1x∗� ,

λ9 ¼ � 1 2D1

0

BBBBBB@

�

þf 0 1 s ∗ ð ÞD1x∗� :

stable.

E2 ζs;

1�ζs D1 ; 0 � � in Rþ

where bx ¼ 1�bs

Lyapunov function L ¼ 1

D1 f 0 1 s

> 0 B@

D1 f 0 1 s

> 0 B@

D1

vuuuuuuuut

D1

vuuuuuuuut

∗ ð Þx∗ þ Q þ R þ V � Z þ 2 � � � f 1 s ∗ ð Þ �

∗ ð Þx∗ þ Q þ R þ V � Z þ 2 � � � f 1 s ∗ ð Þ �

Q2 þ 2QR þ 2QV � 2QZ � 2Q f 1 s

These results lead to the following theorem.

use the similar technique as in [3, 18].

(i) Global stability analysis of E2 ζs;

Consider system (3) restricted to Rþ

nents of E2ð Þ bs; bx , i.e., bs ¼ ζs and bx ¼ 1�ζs

positive definite function with respect to E2 in Rþ

By differentiating L with respect to time t, we get

D1 . From (5), we have

2.4. Global stability and uniform persistence analysis

2 n1ð Þ s �bs

dL

Now we will analyse the global stability of the equilibrium point E2 ζs;

1�ζs D1 ; 0 � �

2 þ n2ð Þ x � bx

dt ¼ n1ð Þ s �bs s

Q2 þ 2QR þ 2QV � 2QZ � 2Q f 1 s

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation…

f 1 s

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

f 1 s

∗ ð Þx∗ 2Q þ 2R þ 2V � 2Z � f

∗ ð Þ 2D1R � 2D1V � 2D1Z � f

∗ ð Þx∗ 2Q þ 2R þ 2V � 2Z � f

∗ ð Þ 2D1R � 2D1V � 2D1Z � f

sx as in (5). Let N be a neighbourhood of equilibrium point

D1 , while n1 and n2 are positive constants. Note that L is a

sx and a Lyapunov function for system (5) in N .

, (7)

sx. To analyse the global stability of the equilibrium point E2, a suitable

0 þ n2ð Þ x � bx x0

2 � � is used, where bs and bx denote the compo-

0 1 s

http://dx.doi.org/10.5772/intechopen.71624

2RV � 2RZ þ V2 � 2VZ þ Z2

0 1 s

2RV � 2RZ þ V2 � 2VZ þ Z2

0 1 s

ð Þ 1�ζs D1 ; 0 � � and the

0 1 s þ

1 CAþ 1

57

CCCCCCA

∗ ð Þ2

þ

1 CAþ 1

CCCCCCA

∗ ð Þ2

∗ ð Þx∗ � � þ R2

∗ ð Þx∗ � � þ f 1 s

∗ ð Þx∗ � � þ R2

∗ ð Þx∗ � � þ f 1 s

∗ ð Þ� f 0 1 s

∗ ð Þ� f 0 1 s

Theorem 7. Suppose that λ8 < 0 and λ9 < 0, then the equilibrium point E3 is locally asymptotically

existence of positive interior equilibrium point E3 s∗; x∗; y∗ ð Þ. For global stability analysis, we

where

$$\begin{aligned} c\_4 &= f\_1'(s^\*) \mathbf{x}^\* + Q + R + V - Z + 3 - \frac{f\_1(s^\*)}{D\_1}, \\ c\_5 &= f\_1'(s^\*) \mathbf{x}^\* (Q + R + V - Z + 2) + \frac{f\_1(s^\*)(Z - R - 2)}{D\_1} + 2Q + 2R + 2V - 2Z + 3, \\ c\_6 &= f\_1'(s^\*) \mathbf{x}^\* (Q + R + V - Z + 1) + \frac{f\_1(s^\*)(Z - R - 1)}{D\_1} + Q + R + V - Z + 1, \end{aligned}$$

and

$$\begin{aligned} Q &= f\_2(\mathbf{x}^\*) y^\* \left[ \left( \frac{D\_1}{1 + \mathbf{x}^\* + D\_1 by^\*} \right) - \left( \frac{D\_1 \mathbf{x}^\*}{(1 + \mathbf{x}^\* + D\_1 by^\*)^2} \right) \right], \\ R &= \left( \frac{f\_2(\mathbf{x}^\*) D\_1^2 b \mathbf{x}^\* y^\*}{(1 + \mathbf{x}^\* + D\_1 by^\*)^2} \right), \\ V &= f\_2'(\mathbf{x}^\*) y^\* \left( \frac{D\_1 \mathbf{x}^\*}{1 + \mathbf{x}^\* + D\_1 by^\*} + 1 \right), \\ Z &= f\_2(\mathbf{x}^\*) \left( \frac{D\_1 \mathbf{x}^\*}{1 + \mathbf{x}^\* + D\_1 by^\*} + 1 \right). \end{aligned}$$

From (6), we obtain the eigenvalues; λ7 ¼ �1, and λ8, λ9 are

Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation… http://dx.doi.org/10.5772/intechopen.71624 57

$$\begin{aligned} \boldsymbol{\theta}\_{8} &= -\frac{1}{2D\_{1}}\Bigl(D\_{1}\big(f\_{1}^{\star}(\mathbf{s}^{\star})\mathbf{x}^{\star} + Q + \mathcal{R} + V - \mathcal{Z} + 2\big) - f\_{1}(\mathbf{s}^{\star})\Bigr) \\ &+ \begin{pmatrix} \left(D\_{1}\left(Q^{2} + 2Q\mathcal{R} + 2QV - 2Q\mathcal{Z} - 2Qf\_{1}(\mathbf{s}^{\star}) - f\_{1}^{\star}(\mathbf{s}^{\star})\mathbf{x}^{\star}\big(2Q + 2\mathcal{R} + 2V - 2\mathcal{Z} - f\_{1}^{\star}(\mathbf{s}^{\star})\mathbf{x}^{\star}\big) + R^{2}\right) + \mathbf{1}\_{\star}\right) \\ &\mathcal{R}V - 2R\mathcal{Z} + V^{2} - 2\mathcal{V}\mathcal{Z} + \mathcal{Z}^{2} \Biggr) \\ &+ f\_{1}(\mathbf{s}^{\star})D\_{1}\mathbf{x}^{\star}\bigr), \\ &\mathcal{I}\_{1}(\mathbf{s}^{\star})D\_{1}\mathbf{x}^{\star}\bigr), \\ \boldsymbol{\theta}\_{8} &= -\frac{1}{2D\_{1}}\Bigl(D\_{1}(f\_{1}^{\star}(\mathbf{s}^{\star})\mathbf{x}^{\star} + Q + \mathcal{R} + V - \mathcal{Z} + 2\right) - f\_{1}(\mathbf{s}^{\star}) \\ &- \begin{pmatrix} \left(D\_{1}\left(\boldsymbol{\big$$

These results lead to the following theorem.

" # !

 Suppose that λ8 < 0 and λ9 < 0, then the equilibrium point E3 is locally asymptotically stable.

#### 2.4. Global stability and uniform persistence analysis

Now we will analyse the global stability of the equilibrium point E2 ζs; ð Þ 1�ζs D1 ; 0 � � and the existence of positive interior equilibrium point E3 s∗; x∗; y∗ ð Þ. For global stability analysis, we use the similar technique as in [3, 18].

(i) Global stability analysis of E2 ζs; 1�ζs D1 ; 0 � �

Consider system (3) restricted to Rþ sx as in (5). Let N be a neighbourhood of equilibrium point E2 ζs; 1�ζs D1 ; 0 � � in Rþ sx. To analyse the global stability of the equilibrium point E2, a suitable Lyapunov function L ¼ 1 2 n1ð Þ s �bs 2 þ n2ð Þ x � bx 2 � � is used, where bs and bx denote the components of E2ð Þ bs; bx , i.e., bs ¼ ζs and bx ¼ 1�ζs D1 , while n1 and n2 are positive constants. Note that L is a positive definite function with respect to E2 in Rþ sx and a Lyapunov function for system (5) in N . By differentiating L with respect to time t, we get

$$\frac{d\mathbf{L}}{dt} = n\_1(\mathbf{s} - \widehat{\mathbf{s}})\mathbf{s}' + n\_2(\mathbf{x} - \widehat{\mathbf{x}})\mathbf{x}',\tag{7}$$

where bx ¼ 1�bs D1 . From (5), we have

$$1 = \widehat{\mathbf{s}} + \widehat{\mathbf{x}} f\_1(\widehat{\mathbf{s}}) \text{ and } D\_1 = f\_1(\widehat{\mathbf{s}}).$$

Hence (7) can be written as

$$\begin{split} L' &= n\_1(\mathbf{s} - \widehat{\mathbf{s}}) \left( \widehat{\mathbf{s}} + \widehat{\mathbf{x}} f\_1(\widehat{\mathbf{s}}) - \mathbf{s} - f\_1(\mathbf{s}) \mathbf{x} \right) + n\_2(\mathbf{x} - \widehat{\mathbf{x}}) \left( \frac{f\_1(\mathbf{s})}{f\_1(\widehat{\mathbf{s}})} - 1 \right) \mathbf{x} \\ &= -n\_1(\mathbf{s} - \widehat{\mathbf{s}})^2 + n\_1(\mathbf{s} - \widehat{\mathbf{s}}) \left( f\_1(\widehat{\mathbf{s}}) \widehat{\mathbf{x}} - f\_1(\mathbf{s}) \mathbf{x} \right) + n\_2(\mathbf{x} - \widehat{\mathbf{x}}) \mathbf{x} \left( \frac{f\_1(\mathbf{s})}{f\_1(\widehat{\mathbf{s}})} - 1 \right) \\ &= n\_{11}(\mathbf{s} - \widehat{\mathbf{s}})^2 + \frac{1}{2} n\_{12}(\mathbf{s} - \widehat{\mathbf{s}}) (\mathbf{x} - \widehat{\mathbf{x}}) + \frac{1}{2} n\_{21}(\mathbf{s} - \widehat{\mathbf{s}}) (\mathbf{x} - \widehat{\mathbf{x}}) + n\_{22}(\mathbf{x} - \widehat{\mathbf{x}})^2. \end{split}$$

where

$$n\_{11} = -n\_1 < 0,$$

$$n\_{12} = n\_{21} = n\_1 \frac{f\_1(\widehat{\mathbf{s}})\widehat{\mathbf{x}} - f\_1(\mathbf{s})\mathbf{x}}{\mathbf{x} - \widehat{\mathbf{x}}},$$

$$n\_{22} = n\_2 \frac{\mathbf{x}(f\_1(\mathbf{s}) - f\_1(\widehat{\mathbf{s}}))}{(\mathbf{x} - \widehat{\mathbf{x}})f\_1(\widehat{\mathbf{s}})}.$$

Clearly that L0 can be written as L0 ¼ XTNX, which T denotes the transpose and the matrix N is particularly a real, symmetric 2 � 2 matrix, where X and N can be represented by

$$X = \begin{pmatrix} \upsilon\_1 \\ \upsilon\_2 \end{pmatrix} = \begin{pmatrix} s - \widehat{s} \\ \varkappa - \widehat{\chi} \end{pmatrix} \text{ and } N = \begin{pmatrix} n\_{11} & \frac{1}{2} n\_{12} \\ 1 & n\_{21} \\ \frac{1}{2} n\_{21} & n\_{22} \end{pmatrix}.$$

Thus, it leads to the following theorem.

Theorem 8. The equilibrium point E2 is global asymptotically stable with respect to solution trajectories are initiated from int Rþ sx if the assumptions n22 < 0 and det N > 0 are satisfied.

Proof. By using the Frobenius Theorem in ([18], Lemma 6.2), we can see that n22 and det N are the leading principal minors of the matrix N: It is shown that matrix N is negative definite if

$$n\_{22} < 0, \text{ and } \det N = \det \begin{pmatrix} n\_{11} & \frac{1}{2} n\_{12} \\ 1 \\ \frac{1}{2} n\_{21} & n\_{22} \end{pmatrix} > 0.1$$

This completes the proof of the theorem. ∎

(ii) Existence of positive interior equilibrium point E3

In this subsection, we present some results of persistence analysis, including uniform persistence and state the necessary conditions for the existence of positive equilibrium point E3. The following lemma from Ref. [19] is applied to obtain the persistence results.

Lemma 2. [19] Let G be an isolated hyperbolic equilibrium point in the omega limit set, ωð Þ X of

Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation…

Proof. This proof strictly depends on Lemma 2. Suppose Η is the region as stated in Theorem 1. It showed that region Η is positive invariant set and any solutions of system (3) emanating at

> ð Þ 1�ζs D1 ; 0

ð Þ E2 \ Ef g2 such that J

þ: Thus, the equilibrium point E2 is not in the omega limit set of P; E2∉ωð Þ P : Next, we

the closure of the orbit of the point J, i.e., clð Þ O Jð Þ must either contains E2 or unbounded. This is

that the boundary flow is isolated and a cyclic with respect to region Η. Thus, the system (3) undergoes uniform persistence and implies that a positive interior equilibrium point

In this section, we investigate Hopf bifurcation on the system (3) with a bifurcation real parameter, σ: Particularly, σ is selected in such a way that the growth rate function f 2 is a

ð Þ E2 ; Ws

is empty which is a contradiction for the positive invariance property of

þ is uniformly bounded. Despites that, the only compact invariant set on ∂R3

þ and J

ð Þ G and Wuð Þ G denote stable and unstable manifolds of G.

is a hyperbolic saddle point and is repels locally in y direction (as in Theorem 5),

is globally asymptotically stable with respect to Rþ

� in ωð Þ X with J

http://dx.doi.org/10.5772/intechopen.71624

þ, i.e., P∈ int R3

¼ ∅, and for unstable mani-

: This completes the proof. ∎

1 ∈ ωð Þ P by Lemma 2. But

þ and J ∈ ωð Þ P : Then

ð Þ 1�ζs D1 ; 0

sx indicates

þ belongs to ωð Þ P , the omega limit set of P: Now,

∉ωð Þ P : Assume that E2 ∈ ωð Þ P is true.

þ ∩ ωð Þ¼ P ∅: We deduce that if E2 ζs;

with respect to Rþ

þ

þ

þ ∩ ωð Þ P 6¼ ∅, and let J ∈ ∂R3

6¼ ∅. Therefore, the result of uniform persistence follows

ð Þ 1�ζs D1 ; 0

ð Þ E2 ∩ int R3

þ

sxy are uniformly bounded and

sx.

þ ∈ Ws

ð Þ G

59

þ is

þ. We shall

the orbit Oð Þ X : Then either ωð Þ¼ X G or there exist points J

ii. system (3) is dissipative and all solutions with initial values in Rþ

. Let P ¼ E3 s∗; x∗; y∗ ð Þ belongs to the interior of R3

þ 1 ∈ Ws

þ ∩ ωð Þ¼ P ∅: Assume that ∂R3

þ

since the omega limit set of P, ωð Þ P must be in int R3

Remark 2. Global stability of equilibrium point E2 ζs;

function of x and σ: Therefore, system (3) takes of the form

þ where ∂R3

ð Þ 1�ζs D1 ; 0

attracted into region Η (as in Theorem 1), and

and J

i. E2 ζs;

a point in R3

ð Þ 1�ζs D1 ; 0

ð Þ E2 ∩ R3

shall show ∂R3

E2 ζs;

Ws

Η ⊂ R3

� ∈ Wuð Þ G , where Ws

ð Þ 1�ζs D1 ; 0

iii. equilibrium point E2 ζs;

Then, the system (3) is uniformly persistence.

show that there is no points Ji ∈∂R3

Then, there is a point J

þ\ Ef g2

fold Wuð Þ E2 ; Wuð Þ E2 ∩ int R3

E3 s∗; x∗; y∗ ð Þ exists (see [20]).

2.5. Hopf bifurcation

we prove that the equilibrium point E2 ζs;

a contradiction, and hence it is proved that ∂R3

is unstable, then, for stable manifold Ws

Theorem 9. Assume that

Lemma 2. [19] Let G be an isolated hyperbolic equilibrium point in the omega limit set, ωð Þ X of the orbit Oð Þ X : Then either ωð Þ¼ X G or there exist points J þ and J � in ωð Þ X with J þ ∈ Ws ð Þ G and J � ∈ Wuð Þ G , where Ws ð Þ G and Wuð Þ G denote stable and unstable manifolds of G.

Theorem 9. Assume that

1 ¼ bs þ bxf 1ð Þ bs and D1 ¼ f 1ð Þ bs :

2 þ n1ð Þ s �bs f 1ð Þ bs bx � f 1ð Þs x � � þ n2ð Þ x � bx x f 1ð Þs

1 2

n11 ¼ �n1 < 0,

n12 ¼ n21 ¼ n1

Clearly that L0 can be written as L0 ¼ XTNX, which T denotes the transpose and the matrix N is

Theorem 8. The equilibrium point E2 is global asymptotically stable with respect to solution trajecto-

Proof. By using the Frobenius Theorem in ([18], Lemma 6.2), we can see that n22 and det N are the leading principal minors of the matrix N: It is shown that matrix N is negative definite if

This completes the proof of the theorem. ∎

In this subsection, we present some results of persistence analysis, including uniform persistence and state the necessary conditions for the existence of positive equilibrium point E3. The

sx if the assumptions n22 < 0 and det N > 0 are satisfied.

n11

1 2

0 B@ 1 2 n12 1

CA > 0:

n21 n22

particularly a real, symmetric 2 � 2 matrix, where X and N can be represented by

¼ s �bs x � bx � �

n22 < 0, and det N ¼ det

following lemma from Ref. [19] is applied to obtain the persistence results.

n22 ¼ n2

and N ¼

f 1ð Þ bs � 1 � �

n21ð Þ s �bs ð Þþ x � bx n22ð Þ x � bx

f 1ð Þ bs bx � f 1ð Þs x x � bx ,

x f 1ð Þ� s f 1ð Þ bs � � ð Þ x � bx f 1ð Þ bs :

n11

1 2

0 B@ 1 2 n12 1 CA:

n21 n22

x

f 1ð Þ bs � 1 � �

> 2 ,

L0 ¼ n1ð Þ s �bs bs þ bxf 1ð Þ� bs s � f 1ð Þs x � � þ n2ð Þ x � bx f 1ð Þs

n12ð Þ s �bs ð Þþ x � bx

Hence (7) can be written as

where

¼ �n1ð Þ s �bs

¼ n11ð Þ s �bs

58 Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

X ¼ v1 v2 � �

(ii) Existence of positive interior equilibrium point E3

Thus, it leads to the following theorem.

ries are initiated from int Rþ

Then, the system (3) is uniformly persistence.

Proof. This proof strictly depends on Lemma 2. Suppose Η is the region as stated in Theorem 1. It showed that region Η is positive invariant set and any solutions of system (3) emanating at a point in R3 þ is uniformly bounded. Despites that, the only compact invariant set on ∂R3 þ is E2 ζs; ð Þ 1�ζs D1 ; 0 . Let P ¼ E3 s∗; x∗; y∗ ð Þ belongs to the interior of R3 þ, i.e., P∈ int R3 þ. We shall show that there is no points Ji ∈∂R3 þ where ∂R3 þ belongs to ωð Þ P , the omega limit set of P: Now, we prove that the equilibrium point E2 ζs; ð Þ 1�ζs D1 ; 0 ∉ωð Þ P : Assume that E2 ∈ ωð Þ P is true.  But Ws ð Þ E2 ∩ R3 þ\ Ef g2 is empty which is a contradiction for the positive invariance property of Η ⊂ R3 þ: Thus, the equilibrium point E2 is not in the omega limit set of P; E2∉ωð Þ P : Next, we shall show ∂R3 þ ∩ ωð Þ¼ P ∅: Assume that ∂R3 þ ∩ ωð Þ P 6¼ ∅, and let J ∈ ∂R3 þ and J ∈ ωð Þ P : Then the closure of the orbit of the point J, i.e., clð Þ O Jð Þ must either contains E2 or unbounded. This is a contradiction, and hence it is proved that ∂R3 þ ∩ ωð Þ¼ P ∅: We deduce that if E2 ζs; ð Þ 1�ζs D1 ; 0 is unstable, then, for stable manifold Ws ð Þ E2 ; Ws ð Þ E2 ∩ int R3 þ ¼ ∅, and for unstable manifold Wuð Þ E2 ; Wuð Þ E2 ∩ int R3 þ 6¼ ∅. Therefore, the result of uniform persistence follows since the omega limit set of P, ωð Þ P must be in int R3 þ : This completes the proof. ∎

Remark 2. Global stability of equilibrium point E2 ζs; ð Þ 1�ζs D1 ; 0 with respect to Rþ sx indicates that the boundary flow is isolated and a cyclic with respect to region Η. Thus, the system (3) undergoes uniform persistence and implies that a positive interior equilibrium point E3 s∗; x∗; y∗ ð Þ exists (see [20]).

#### 2.5. Hopf bifurcation

In this section, we investigate Hopf bifurcation on the system (3) with a bifurcation real parameter, σ: Particularly, σ is selected in such a way that the growth rate function f 2 is a function of x and σ: Therefore, system (3) takes of the form

$$\begin{cases} \mathbf{s}' = \mathbf{1} - \mathbf{s} - f\_1(\mathbf{s})\mathbf{x}, \\\\ \mathbf{x}' = \left(\frac{f\_1(\mathbf{s})}{D\_1} - \mathbf{1}\right) \mathbf{x} - f\_2(\mathbf{x}; \sigma) y \left(\mathbf{1} + \frac{D\_1 \mathbf{x}}{1 + \mathbf{x} + D\_1 \mathbf{b} y}\right), \\\\ \mathbf{y}' = f\_2(\mathbf{x}; \sigma) y \left(\mathbf{1} + \frac{D\_1 \mathbf{x}}{1 + \mathbf{x} + D\_1 \mathbf{b} y}\right) - y, \end{cases} \tag{8}$$

where sð Þ0 > 0, xð Þ0 > 0, yð Þ0 > 0: Next, we do linearization on the system (8). First, let

$$\begin{cases} \mathbf{S} = \mathbf{s} - h\_1; \\\\ \mathbf{X} = \mathbf{x} - h\_2; \quad \Rightarrow \\\\ \mathbf{Y} = \mathbf{y} - h\_3; \end{cases} \Rightarrow \begin{cases} \mathbf{s} = \mathbf{S} + h\_1; \\\\ \mathbf{x} = \mathbf{X} + h\_2; \\\\ \mathbf{y} = \mathbf{Y} + h\_3; \end{cases}$$

where ð Þ h1; h2; h3 is the non-trivial equilibrium point. Then, we obtain the following differential equations

$$\begin{cases} \begin{aligned} S' &= 1 - S - h\_1 - f\_1(S + h\_1)(X + h\_2), \\ X' &= \left(\frac{f\_1(S + h\_1)}{D\_1} - 1\right)(X + h\_2) - f\_2(X + h\_2; \sigma)(Y + h\_3) \left(1 + \frac{D\_1(X + h\_2)}{1 + (X + h\_2) + D\_1b(Y + h\_3)}\right), \\ &\quad Y' = f\_2(X + h\_2; \sigma)(Y + h\_3) \left(1 + \frac{D\_1(X + h\_2)}{1 + (X + h\_2) + D\_1b(Y + h\_3)}\right) - (Y + h\_3). \end{aligned} \end{cases} \tag{9}$$

The Jacobian matrix of system (8) is given by

$$J\_{4} = \begin{bmatrix} -f\_{1}'(\mathbf{s})\mathbf{x} - \mathbf{1} & -f\_{1}(\mathbf{s})\\ \frac{f\_{1}'(\mathbf{s})\mathbf{x}}{D\_{1}} & \frac{f\_{1}(\mathbf{s})}{D\_{1}} - f\_{2}(\mathbf{x};\sigma)y\left[\left(\frac{D\_{1}}{\mathbf{1}+\mathbf{x}+D\_{1}\mathbf{b}}\right) - \left(\frac{D\_{1}\mathbf{x}}{\left(\mathbf{1}+\mathbf{x}+D\_{1}\mathbf{b}\right)^{2}}\right)\right] - f\_{2}'(\mathbf{x};\sigma)y\left(\left(\frac{D\_{1}\mathbf{x}}{\mathbf{1}+\mathbf{x}+D\_{1}\mathbf{b}}\right) + 1\right) - 1\\ 0 & f\_{2}(\mathbf{x};\sigma)y\left[\left(\frac{D\_{1}}{\mathbf{1}+\mathbf{x}+D\_{1}\mathbf{b}}\mathbf{y}\right) - \left(\frac{D\_{1}\mathbf{x}}{\left(\mathbf{1}+\mathbf{x}+D\_{1}\mathbf{b}\right)^{2}}\right)\right] + f\_{2}'(\mathbf{x};\sigma)y\left(\mathbf{1} + \frac{D\_{1}\mathbf{x}}{\mathbf{1}+\mathbf{x}+D\_{1}\mathbf{b}}\right) \end{bmatrix}$$

$$\begin{aligned} f\_2(\mathbf{x}; \sigma) \left[ \left( \frac{D\_1^2 b \mathbf{x} \mathbf{y}}{\left(1 + \mathbf{x} + D\_1 b \mathbf{y}\right)^2} \right) - \left( \frac{D\_1 \mathbf{x}}{1 + \mathbf{x} + D\_1 b \mathbf{y}} + 1 \right) \right] \\ f\_2(\mathbf{x}; \sigma) \left[ \left( 1 + \frac{D\_1 \mathbf{x}}{1 + \mathbf{x} + D\_1 b \mathbf{y}} \right) - \left( \frac{D\_1^2 b \mathbf{x} \mathbf{y}}{\left(1 + \mathbf{x} + D\_1 b \mathbf{y}\right)^2} \right) \right] - 1 \end{aligned}$$

:

Jσð Þ E2 ¼

where

� f 0

0

BBBBBBBBBBBBBBBB@

f 0

c1 ¼ 1

þ D1f2

c2 ¼ 1

� f 2

� f 2

obtain the matrices

1ð Þ ζs ð Þ 1 � ζs D1

1ð Þ ζs ð Þ 1 � ζs D2 1

The characteristic equation of (10) is given as

D1ð Þ D1 � ζs þ 1 f

1 � ζs D1

D1ð Þ D1 � ζs þ 1 f

1 � ζs D1 ; σ

1 � ζs D1 ; σ � � 2D1 þ 4D2

; σ � �

� f 1ð Þ ζs 2 þ 2D1 � 2ζs þ f 2

c3 ¼ 1

�f 2

�f 2

� 1 �f 1ð Þ ζs 0

Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation…

1 � ζs D1 ; σ � � 1 þ

1 � ζs D1 ; σ � � 1 þ 1 � ζs

http://dx.doi.org/10.5772/intechopen.71624

1 � ζs D1 � �

1

CCA

1

1

1

61

CCCCCCCCCCCCCCCCA :

(10)

CCA

CCA

1

CCA � 1

1 þ

1 � ζs

1 � ζs D1 � �

0

0

BB@

0

BB@

s � D1ζs � � þ f1 ζs ð Þð Þ ζs � D1 � 1 �

ðζs þ D1ζs � 2D1 � 1Þ þ 3D1ð ÞÞ 1 þ D1 � ζs ,

s þ 2 þ 2D1 � 4ζs � 2D1ζs

� ��

� �

1 � ζs D1 ; σ � �

1 � 2D1ζs � 2D2

0

; σ

We applied Routh Hurwitz criterion (see [21, 22]) onto the characteristic equation (11) and

� �

s � 3D1ζs þ D1ζ2

1ζs � � þ 3D1ζs � 3D2

s � D1 þ D1ζ2

1 ζs ð Þ 1 � 2ζs þ D1 þ ζ2

� ��

� � ζs � 1 � 2D1 þ D1ζs ð Þ þ 1 � ζs � �

ð Þþ D1 � D1ζs D1ζs � D1 � D2

BB@

0

BB@

λ3 þ c1λ2 þ c2λ þ c3 ¼ 0, (11)

s

ð Þ ζs � 2D1 � 1 þ D1ζs

s

s � D1ζs

1Þ:

1 � 3D1Þ,

1 þ

� 1 �f 2

1 ζs ð Þ 1 þ D1 � 2ζs þ ζ2

f 1ð Þ ζs D1

0 0 f 2

0

0 1ð Þ ζs 2ζ2

� �

� � 1 þ 2D1 � 2ζs þ ζ2

D1ð Þ D1 � ζs þ 1 f

; σ � � 1 � 2ζs þ ζ2

; σ � �

1 � ζs D1

1 � ζs D1

1 � ζs D1

�f 1 ζs ð Þ f 2

Thus, the Jacobian matrix of system (8) about E2 ζs; ð Þ 1�ζs D1 ; 0; σ � � is Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation… http://dx.doi.org/10.5772/intechopen.71624 61

$$f\_{o(E\_{2})} = \begin{pmatrix} -\frac{f\_{1}'(\zeta\_{s})(1-\zeta\_{s})}{D\_{1}} - 1 & -f\_{1}(\zeta\_{s}) & 0\\ \frac{f\_{1}'(\zeta\_{s})(1-\zeta\_{s})}{D\_{1}^{2}} & \frac{f\_{1}(\zeta\_{s})}{D\_{1}} - 1 & -f\_{2}\left(\frac{1-\zeta\_{s}}{D\_{1}};\sigma\right)\left(1 + \left(\frac{1-\zeta\_{s}}{1+\left(\frac{1-\zeta\_{s}}{D\_{1}}\right)}\right)\right) \\\\ 0 & 0 & f\_{2}\left(\frac{1-\zeta\_{s}}{D\_{1}};\sigma\right)\left(1 + \left(\frac{1-\zeta\_{s}}{1+\left(\frac{1-\zeta\_{s}}{D\_{1}}\right)}\right)\right) - 1\\\\ \end{pmatrix}.\tag{10}$$

The characteristic equation of (10) is given as

$$
\lambda^3 + c\_1 \lambda^2 + c\_2 \lambda + c\_3 = 0,\tag{11}
$$

where

s0 ¼ 1 � s � f 1ð Þs x,

where sð Þ0 > 0, xð Þ0 > 0, yð Þ0 > 0: Next, we do linearization on the system (8). First, let

)

where ð Þ h1; h2; h3 is the non-trivial equilibrium point.

 Suppose that the equilibrium point E2 ζs;

f 1 s∗ ð Þ D1

0 f 2 x∗ ð Þ ; σ y∗ D1

Re dλj dσ � �

1

where

Since K σ∗

and λ3 σ∗ 1 � � ¼ �c1 σ∗

from E2 ζs;

Jσð Þ E3 ¼

ð Þ 1�ζs D1 ; 0; σ

� f 0

> f 0 1 s∗ ð Þx∗ D1

1 � �M σ∗ 1 � � þ L σ∗

and calculating the implicit derivative, we obtain the following equations

ð Þ σ c 0 1ð Þþ σ c 0

� � 6¼ 0, we have

� � for some values of σ in the neighbourhood of σ∗

1 s∗ ð Þx∗ � 1 �f 1 s∗ ð Þ

1 þ x∗ þ D1by∗ � �

1 þ x∗ þ D1by∗ � �

matrix of the system (8) about the equilibrium point E3 is given by

� f 2 x∗ ð Þ ; σ y∗ D1

Hence, the characteristic equation for the Jacobian matrix Jσð Þ E3 is

Nð Þ¼ σ 2α σð Þβ σð Þc

¼ K σ∗ 1 � �M σ∗ 1 � � þ L σ∗

σ¼σ∗ 1

ð Þ� σ Lð Þ σ β0

ð Þþ σ Kð Þ σ β0

Lð Þ¼ σ 6α σð Þβ σð Þþ 2c1ð Þ σ β σð Þ;

K2 σ∗ 1 � � þ L2 σ∗

0 1ð Þþ σ c 0 2ð Þ σ β σð Þ:

By substituting λ1ð Þ¼ σ α σð Þþ β σð Þi and λ2ð Þ¼ σ α σð Þ� β σð Þi into characteristic equation (11)

Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation…

ð Þþ σ Mð Þ¼ σ 0,

ð Þþ σ Nð Þ¼ σ 0,

ð Þþ σ 2c1ð Þ σ α σð Þþ c2ð Þ� σ 3β2

2ð Þ σ α σð Þ� c0

1 � �N σ∗ 1 � �

� � 6¼ 0: We conclude the details above in the following theorem.

similar as in (2) for the system (3) together with hypothesis H1 until H3 hold. Then the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating

Next, we determine the Hopf bifurcation at the equilibrium point E3 s∗; x∗; y∗ ð Þ ; σ : The Jacobian

� D1x∗ 1 þ x∗ þ D1by∗ ð Þ2

� D1x∗ 1 þ x∗ þ D1by∗ ð Þ2

f 2 x∗ ð Þ ; σ D2

f 2 x∗ ð Þ ; σ D1x∗

λ3 þ c4λ2 þ c5λ þ c6 ¼ 0, (16)

" # !

" # !

1

ð Þ 1�ζs D1 ; 0; σ

6¼ 0, j ¼ 1, 2, 3: (14)

http://dx.doi.org/10.5772/intechopen.71624

(15)

63

�

1ð Þ σ β2

� � 6¼ 0, j ¼ 1, 2, 3,

1:

� f 2

þ f 0

1bx∗y∗ 1 þ x∗ þ D1by∗ ð Þ2 � D1x∗

1 þ x∗ þ D1by∗ � D2

0 x∗ ð Þ ; σ y∗ D1x∗

2 x∗ ð Þ ; σ y∗ D1x∗

0

!

!

1 þ x∗ þ D1by∗ þ 1 � �

1 þ x∗ þ D1by∗ þ 1 � �

1 þ x∗ þ D1by∗ � 1

1bx∗y∗ 1 þ x∗ þ D1by∗ ð Þ2 þ 1 � 1

� 1

ð Þ σ ;

ð Þþ σ c0

3ð Þ σ ;

� � exists and those assumptions

When the assumptions on functional responses f 1ð Þ¼ s f 1ð Þ ζs and f 2ð Þ¼ x f 2 1�ζs D1 ; σ � � as in (2) for the system (3), together with the hypotheses H1 until H3;

H1 : D1 � ζs þ 1 > 0, H2 : f 0 1ð Þ ζs 1 þ D1 � 2ζs þ ζ2 s � D1ζs � � þ f 1ð Þ ζs ð Þ ζs � D1 � 1 þ D1f 2 1 � ζs D1 ; σ � �ðζs þ D1ζs � 2D1 � 1Þ þ 3D1ð Þ 1 þ D1 � ζs > 0, H3 : f 0 1ð Þ ζs 1 � 2ζs þ D1 þ ζ2 s � D1ζs � f 2 1 � ζs D1 ; σ � � 1 � 2ζs þ ζ2 s � D1 þ D1ζ2 s � � � � � f 1ð Þ ζs f 2 1 � ζs D1 ; σ � �ðζs � 1 � 2D1 þ D1ζsÞ þ 1 � ζs � � � f 2 1 � ζs D1 ; σ � �ð Þþ D1 � D1ζs D1ζs � D1 � D2 1 > 0,

hold, we will have c1 > 0 and c3 > 0: We shall obtain two pure imaginary roots for the characteristic equation (11) if and only if c1c2 ¼ c3 for some values of σ, say, σ∗ 1:

Since at σ ¼ σ∗ 1, there exists an interval σ∗ 1 � ε; σ∗ 1 þ ε � � containing σ∗ 1 for some ε > 0 such that σ∈ σ∗ 1 � ε; σ∗ 1 þ ε � �. Thus, for σ∈ σ∗ 1 � ε; σ∗ 1 þ ε � �, the characteristic equation (11) cannot have positive real roots. For σ ¼ σ∗ 1, we acquire (see [3, 10])

$$(\lambda^2 + c\_2)(\lambda + c\_1) = 0,\tag{13}$$

that consist of three roots; λ1 ¼ ffiffiffiffi c2 p i, λ2 ¼ � ffiffiffiffi c2 p i, and λ3 ¼ �c1: As for σ∈ σ∗ 1 � ε; σ∗ 1 þ ε � �, all the roots are in general of the form;

$$
\lambda\_1(\sigma) = \alpha(\sigma) + \beta(\sigma)i,
$$

$$
\lambda\_2(\sigma) = \alpha(\sigma) - \beta(\sigma)i,
$$

$$
\lambda\_3(\sigma) = -\mathbf{c}\_1(\sigma).
$$

In order to apply the Hopf bifurcation theorem as stated in [9, 23] towards the system (8), we must verify the transversality condition

Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation… http://dx.doi.org/10.5772/intechopen.71624 63

$$\operatorname{Re}\left[\frac{d\lambda\_j}{d\sigma}\right]\_{\sigma=\sigma\_1^\*} \neq 0, \qquad j=1,2,3. \tag{14}$$

By substituting λ1ð Þ¼ σ α σð Þþ β σð Þi and λ2ð Þ¼ σ α σð Þ� β σð Þi into characteristic equation (11) and calculating the implicit derivative, we obtain the following equations

$$\begin{aligned} K(\sigma)\alpha'(\sigma) - L(\sigma)\beta'(\sigma) + M(\sigma) &= 0, \\ L(\sigma)\alpha'(\sigma) + K(\sigma)\beta'(\sigma) + N(\sigma) &= 0, \end{aligned} \tag{15}$$

where

M1 ¼ ½ � c1 ; M2 ¼ c1 1

62 Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

c3 c2 � �

Thus, the characteristic equation (11) has all negative real parts of λ if and only if

When the assumptions on functional responses f 1ð Þ¼ s f 1ð Þ ζs and f 2ð Þ¼ x f 2

for the system (3), together with the hypotheses H1 until H3;

1 � ζs D1 ; σ � �

1 � ζs D1 ; σ � �

1, there exists an interval σ∗

s � D1ζs � � þ f 1ð Þ ζs ð Þ ζs � D1 � 1

s � D1ζs � f 2

� �

ð Þþ D1 � D1ζs D1ζs � D1 � D2

characteristic equation (11) if and only if c1c2 ¼ c3 for some values of σ, say, σ∗

1 � ε; σ∗

1, we acquire (see [3, 10])

λ2 þ c2

p i, λ2 ¼ � ffiffiffiffi

c2

H1 : D1 � ζs þ 1 > 0,

1ð Þ ζs 1 þ D1 � 2ζs þ ζ2

1ð Þ ζs 1 � 2ζs þ D1 þ ζ2

1 � ζs D1 ; σ � �

1 þ ε � �.  For σ ¼ σ∗

that consist of three roots; λ1 ¼ ffiffiffiffi

the roots are in general of the form;

must verify the transversality condition

σ∈ σ∗

; M3 ¼

c1 > 0 c3 > 0 9 >=

>;

ðζs þ D1ζs � 2D1 � 1Þ þ 3D1ð Þ 1 þ D1 � ζs > 0,

1 � 2ζs þ ζ2

1 þ ε � �, the characteristic equation (11) cannot have

� �ð Þ¼ λ þ c1 0, (13)

p i, and λ3 ¼ �c1: As for σ∈ σ∗

1 > 0,

1 þ ε � � containing σ∗

c1c2 � c3 > 0:

1 � ζs D1 ; σ � �

� � � �

ðζs � 1 � 2D1 þ D1ζsÞ þ 1 � ζs

hold, we will have c1 > 0 and c3 > 0: We shall obtain two pure imaginary roots for the

1 � ε; σ∗

c2

λ1ð Þ¼ σ α σð Þþ β σð Þi,

λ2ð Þ¼ σ α σð Þ� β σð Þi,

λ3ð Þ¼� σ c1ð Þ σ :

In order to apply the Hopf bifurcation theorem as stated in [9, 23] towards the system (8), we

2 6 4

c1 1 0 c3 c2 c1 0 0 c3 3 7 5:

(12)

as in (2)

1�ζs D1 ; σ � �

s � D1 þ D1ζ2

s

1:

1 for some ε > 0 such that

1 � ε; σ∗ 1 þ ε � �, all

$$\begin{aligned} K(\sigma) &= 3\alpha^2(\sigma) + 2c\_1(\sigma)\alpha(\sigma) + c\_2(\sigma) - 3\beta^2(\sigma); \\\\ M(\sigma) &= \alpha^2(\sigma)c\_1'(\sigma) + c\_2'(\sigma)\alpha(\sigma) - c\_1'(\sigma)\beta^2(\sigma) + c\_3'(\sigma); \\\\ L(\sigma) &= 6\alpha(\sigma)\beta(\sigma) + 2c\_1(\sigma)\beta(\sigma); \\\\ N(\sigma) &= 2\alpha(\sigma)\beta(\sigma)c\_1'(\sigma) + c\_2'(\sigma)\beta(\sigma). \end{aligned}$$

Since K σ∗ 1 � �M σ∗ 1 � � þ L σ∗ 1 � �N σ∗ 1 � � 6¼ 0, we have

$$\operatorname{Re}\left[\frac{d\lambda\_{\rangle}}{d\sigma}\right]\_{\sigma=\sigma\_1^\*} = \frac{K\left(\sigma\_1^\*\right)M\left(\sigma\_1^\*\right) + L\left(\sigma\_1^\*\right)N\left(\sigma\_1^\*\right)}{K^2\left(\sigma\_1^\*\right) + L^2\left(\sigma\_1^\*\right)} \neq 0, \quad j = 1, 2, 3\ldots$$

and λ3 σ∗ 1 � � ¼ �c1 σ∗ 1 � � 6¼ 0: We conclude the details above in the following theorem. Theorem 10. Suppose that the equilibrium point E2 ζs; ð Þ 1�ζs D1 ; 0; σ � � exists and those assumptions similar as in (2) for the system (3) together with hypothesis H1 until H3 hold. Then the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating

$$\text{from } E\_2\left(\zeta\_s, \frac{(1-\zeta\_s)}{D\_1}, 0; \sigma\right) \text{ for some values of } \sigma \text{ in the neighborhood of } \sigma\_1^\*.$$

Next, we determine the Hopf bifurcation at the equilibrium point E3 s∗; x∗; y∗ ð Þ ; σ : The Jacobian matrix of the system (8) about the equilibrium point E3 is given by

Jσð Þ E3 ¼ � f 0 1 s∗ ð Þx∗ � 1 �f 1 s∗ ð Þ f 0 1 s∗ ð Þx∗ D1 f 1 s∗ ð Þ D1 � f 2 x∗ ð Þ ; σ y∗ D1 1 þ x∗ þ D1by∗ � � � D1x∗ 1 þ x∗ þ D1by∗ ð Þ2 " # ! � f 2 0 x∗ ð Þ ; σ y∗ D1x∗ 1 þ x∗ þ D1by∗ þ 1 � � � 1 0 f 2 x∗ ð Þ ; σ y∗ D1 1 þ x∗ þ D1by∗ � � � D1x∗ 1 þ x∗ þ D1by∗ ð Þ2 " # ! þ f 0 2 x∗ ð Þ ; σ y∗ D1x∗ 1 þ x∗ þ D1by∗ þ 1 � � 2 6 6 6 6 6 6 6 4 0 f 2 x∗ ð Þ ; σ D2 1bx∗y∗ 1 þ x∗ þ D1by∗ ð Þ2 � D1x∗ 1 þ x∗ þ D1by∗ � 1 ! f 2 x∗ ð Þ ; σ D1x∗ 1 þ x∗ þ D1by∗ � D2 1bx∗y∗ 1 þ x∗ þ D1by∗ ð Þ2 þ 1 ! � 1 3 7 7 7 7 7 7 7 5 :

Hence, the characteristic equation for the Jacobian matrix Jσð Þ E3 is

$$
\lambda^3 + c\_4 \lambda^2 + c\_5 \lambda + c\_6 = 0,\tag{16}
$$

where

$$\begin{aligned} c\_4 &= f\_1'(s^\*) \mathbf{x}^\* + Q + R + V - Z + 3 - \frac{f\_1(s^\*)}{D\_1}, \\ c\_5 &= f\_1'(s^\*) \mathbf{x}^\* (Q + R + V - Z + 2) + \frac{f\_1(s^\*)(Z - R - 2)}{D\_1} + 2Q + 2R + 2V - 2Z + 3, \\ c\_6 &= f\_1'(s^\*) \mathbf{x}^\* (Q + R + V - Z + 1) + \frac{f\_1(s^\*)(Z - R - 1)}{D\_1} + Q + R + V - Z + 1, \end{aligned}$$

and

$$\begin{aligned} Q &= f\_2(\mathbf{x}^\*; \sigma) y^\* \left[ \left( \frac{D\_1}{1 + \mathbf{x}^\* + D\_1 b y^\*} \right) - \left( \frac{D\_1 \mathbf{x}^\*}{(1 + \mathbf{x}^\* + D\_1 b y^\*)^2} \right) \right], \\ R &= \left( \frac{f\_2(\mathbf{x}^\*; \sigma) D\_1^2 b \mathbf{x}^\* y^\*}{(1 + \mathbf{x}^\* + D\_1 b y^\*)^2} \right), \\ V &= f\_2'(\mathbf{x}^\*; \sigma) y^\* \left( \frac{D\_1 \mathbf{x}^\*}{1 + \mathbf{x}^\* + D\_1 b y^\*} + 1 \right), \\ Z &= f\_2(\mathbf{x}^\*; \sigma) \left( \frac{D\_1 \mathbf{x}^\*}{1 + \mathbf{x}^\* + D\_1 b y^\*} + 1 \right). \end{aligned}$$

By applying the Routh-Hurwitz criterion (see [21, 22]) towards the characteristic equation (16), we obtain the following Hurwitz matrices

$$M\_4 = [c\_4];\ M\_5 = \begin{bmatrix} c\_4 & 1\\ c\_6 & c\_5 \end{bmatrix};\ M\_6 = \begin{bmatrix} c\_4 & 1 & 0\\ c\_6 & c\_5 & c\_4\\ 0 & 0 & c\_6 \end{bmatrix}.$$

Thus, the characteristic equation (16) has all negative real parts of λ if and only if

$$\begin{cases} c\_4 > 0 \\ c\_6 > 0 \\ c\_4c\_5 - c\_6 > 0. \end{cases} \tag{17}$$

λ2 þ c5

p i, λ2 ¼ � ffiffiffiffi

σ¼σ∗ 2

ð Þ� σ A2ð Þ σ β0

ð Þþ σ A1ð Þ σ β0

0 1ð Þþ σ c 0 2ð Þ σ β σð Þ:

� � 6¼ 0, we have

∗; 1 � s∗

exists and those assumptions similar as in (2) for the system (3) together with hypotheses H4 until H6 hold. Then the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating from E3 s∗; x∗; y∗ ð Þ ; σ for some values of σ in the neighbourhood of σ∗

We have proposed and analysed a simple nutrient-predator-prey model in a chemostat with intratrophic predation. This system consisted of the nutrient s, prey organisms x and predator organisms y. We conclude that intratrophic predation denoted as b does not affected the local stability and existence of E1 and E2 when no predator organisms involved. Next, we have

A1 σ∗ 2 � � � � 2 þ A2 σ∗

c5

Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation…

λ1ð Þ¼ σ α σð Þþ β σð Þi, λ2ð Þ¼ σ α σð Þ� β σð Þi, λ3ð Þ¼� σ c4ð Þ σ ,

By substituting λ1ð Þ¼ σ α σð Þþ β σð Þi and λ2ð Þ¼ σ α σð Þ� β σð Þi into characteristic equation (18)

ð Þþ σ 2c1ð Þ σ α σð Þþ c2ð Þ� σ 3β2

2ð Þ σ α σð Þ� c

2 � �B2 σ∗ 2 � �

2

� � 6¼ 0: We summarise the discussion above in the following theorem.

f 1 s∗ ð Þ ; 2x ð Þ ∗ þ 1 f 2 x∗ ð Þ� x∗ � 1

!

ð Þþ σ B1ð Þ¼ σ 0,

ð Þþ σ B2ð Þ¼ σ 0,

0 1ð Þ σ β2

c5

To establish Hopf bifurcation towards system (8), we must show that

and calculating the implicit derivative, we get the following equations;

A2ð Þ¼ σ 6α σð Þβ σð Þþ 2c1ð Þ σ β σð Þ;

ð Þ σ c 0 1ð Þþ σ c 0

B2ð Þ¼ σ 2α σð Þβ σð Þc

¼ A1 σ∗ 2 � �B1 σ∗ 2 � � þ A2 σ∗

∗; x∗; y∗ ð Þ¼ ; σ E3 s

Re dλj dσ � �

A1ð Þ σ α0

A2ð Þ σ α0

A1ð Þ¼ σ 3α2

B1ð Þ¼ σ α2

2 � �B2 σ∗ 2

σ¼σ∗ 2

Consisting of three roots; λ1 ¼ ffiffiffiffi

roots are in general of the form;

where

Since A1 σ∗

and λ3 σ∗ 2 � � ¼ �c4 σ∗

2.6. Discussion

2 � �B1 σ∗ 2 � � þ A2 σ∗

> Re dλj dσ � �

> > 2

E3 s

Theorem 11. Suppose that the equilibrium point

� �ð Þ¼ λ þ c4 0, (18)

6¼ 0, j ¼ 1, 2, 3: (19)

�

ð Þ σ ;

ð Þþ σ c 0 3ð Þ σ ;

� � � � 2 6¼ 0, j ¼ 1, 2, 3,

bD1 1 � f 2 x∗ ð Þ � � ; σ

2 � ε; σ∗ 2 þ ε � �, all

http://dx.doi.org/10.5772/intechopen.71624

(20)

65

2:

p i, λ3 ¼ �c4: As for σ∈ σ∗

Suppose the assumptions of functional response f 1 s∗ ð Þ and f 2 x∗ ð Þ ; σ similar as in (2) for the system (3), together with the hypotheses H4 until H6;

$$\begin{aligned} H\_4: &f\_1'(s^\*) \mathbf{x}^\* + Q + \mathbf{R} + V + \mathbf{3} > \frac{f\_1(s^\*)}{D\_1} + \mathbf{Z}, \\ H\_5: &Q + \mathbf{R} + V + \mathbf{1} > \mathbf{Z}, \\ H\_6: &\mathbf{Z} > \mathbf{R} + \mathbf{1}, \end{aligned}$$

hold, then clearly c4 > 0 and c6 > 0. In particular, we shall have two pure imaginary roots for the characteristic equation (16) if and only if c4c5 ¼ c6 for some values of σ, say, σ∗ 2: Since at σ ¼ σ∗ 2, there exists an interval σ∗ 2 � ε; σ∗ 2 þ ε � � containing σ∗ 2 for some ε > 0. Then, for σ∈ σ∗ 2 � ε; σ∗ 2 � þεÞ, the characteristic Eq. (16) cannot have positive real roots. For σ ¼ σ∗ 2, we get (see [3, 10])

Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation… http://dx.doi.org/10.5772/intechopen.71624 65

$$(\lambda^2 + c\_5)(\lambda + c\_4) = 0,\tag{18}$$

Consisting of three roots; λ1 ¼ ffiffiffiffi c5 p i, λ2 ¼ � ffiffiffiffi c5 p i, λ3 ¼ �c4: As for σ∈ σ∗ 2 � ε; σ∗ 2 þ ε � �, all roots are in general of the form;

$$\begin{array}{l} \lambda\_1(\sigma) = \alpha(\sigma) + \beta(\sigma)\dot{\iota}, \\ \lambda\_2(\sigma) = \alpha(\sigma) - \beta(\sigma)\dot{\iota}, \\ \lambda\_3(\sigma) = -c\_4(\sigma), \end{array}$$

To establish Hopf bifurcation towards system (8), we must show that

$$\operatorname{Re}\left[\frac{d\lambda\_j}{d\sigma}\right]\_{\sigma=\sigma\_2^\*} \neq 0, \qquad j=1,2,3. \tag{19}$$

By substituting λ1ð Þ¼ σ α σð Þþ β σð Þi and λ2ð Þ¼ σ α σð Þ� β σð Þi into characteristic equation (18) and calculating the implicit derivative, we get the following equations;

$$\begin{cases} A\_1(\sigma)\alpha'(\sigma) - A\_2(\sigma)\beta'(\sigma) + B\_1(\sigma) = 0, \\ A\_2(\sigma)\alpha'(\sigma) + A\_1(\sigma)\beta'(\sigma) + B\_2(\sigma) = 0, \end{cases} \tag{20}$$

where

where

and

c4 ¼ f 0 1 s

c5 ¼ f 0 1 s

c6 ¼ f 0 1 s ∗ ð Þx∗

∗ ð Þx∗ þ Q þ R þ V � Z þ 3 � f 1 s∗ ð Þ

64 Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

Q ¼ f 2 x∗ ð Þ ; σ y∗ D1

1 þ x∗ þ D1by∗ ð Þ2 !

2 x∗ ð Þ ; σ y∗ D1x∗

M4 ¼ ½ � c4 ; M5 ¼ c4 1

R ¼ f 2 x∗ ð Þ ; σ D2

V ¼ f 0

we obtain the following Hurwitz matrices

Z ¼ f 2 x∗ ð Þ ; σ

system (3), together with the hypotheses H4 until H6;

there exists an interval σ∗

H4 : f 0 1 s

2 � ε; σ∗

∗ ð Þx∗ðQ þ R þ V � Z þ 2Þ þ f 1 s∗ ð Þð Þ Z � R � 2

ðQ þ R þ V � Z þ 1Þ þ f 1 s∗ ð Þð Þ Z � R � 1

1 þ x∗ þ D1by∗ � �

,

1 þ x∗ þ D1by∗ þ 1 � �

By applying the Routh-Hurwitz criterion (see [21, 22]) towards the characteristic equation (16),

c6 c5 � �

Thus, the characteristic equation (16) has all negative real parts of λ if and only if

H5 : Q þ R þ V þ 1 > Z,

characteristic equation (16) if and only if c4c5 ¼ c6 for some values of σ, say, σ∗

2 þ ε � � containing σ∗

H6 : Z > R þ 1,

þεÞ, the characteristic Eq. (16) cannot have positive real roots. For σ ¼ σ∗

1bx∗y∗

D1x∗ 1 þ x∗ þ D1by∗ þ 1 � �

D1 ,

D1

D1

" # !

,

2 6 4

9 >=

>;

:

; M6 ¼

c4 > 0 c6 > 0 c4c5 � c6 > 0:

∗ ð Þx∗ þ Q þ R þ V þ 3 > f 1 s∗ ð Þ

Suppose the assumptions of functional response f 1 s∗ ð Þ and f 2 x∗ ð Þ ; σ similar as in (2) for the

hold, then clearly c4 > 0 and c6 > 0. In particular, we shall have two pure imaginary roots for the

� D1x∗

c4 1 0 c6 c5 c4 0 0 c6

D1

þ Z,

2 for some ε > 0. Then, for σ∈ σ∗

3 7 5:

1 þ x∗ þ D1by∗ ð Þ2

þ 2Q þ 2R þ 2V � 2Z þ 3,

,

(17)

2,

2 � ε; σ∗ 2

2: Since at σ ¼ σ∗

�

2, we get (see [3, 10])

þ Q þ R þ V � Z þ 1,

$$\begin{aligned} A\_1(\sigma) &= 3\alpha^2(\sigma) + 2c\_1(\sigma)\alpha(\sigma) + c\_2(\sigma) - 3\beta^2(\sigma); \\ A\_2(\sigma) &= 6\alpha(\sigma)\beta(\sigma) + 2c\_1(\sigma)\beta(\sigma); \\\\ B\_1(\sigma) &= \alpha^2(\sigma)c\_1'(\sigma) + c\_2'(\sigma)\alpha(\sigma) - c\_1'(\sigma)\beta^2(\sigma) + c\_3'(\sigma); \\ B\_2(\sigma) &= 2\alpha(\sigma)\beta(\sigma)c\_1'(\sigma) + c\_2'(\sigma)\beta(\sigma). \end{aligned}$$

Since A1 σ∗ 2 � �B1 σ∗ 2 � � þ A2 σ∗ 2 � �B2 σ∗ 2 � � 6¼ 0, we have

$$\operatorname{Re}\left[\frac{d\lambda\_j}{d\sigma}\right]\_{\sigma=\sigma\_2^\*} = \frac{A\_1(\sigma\_2^\*)B\_1(\sigma\_2^\*) + A\_2(\sigma\_2^\*)B\_2(\sigma\_2^\*)}{\left(A\_1(\sigma\_2^\*)\right)^2 + \left(A\_2(\sigma\_2^\*)\right)^2} \neq 0, \quad j = 1, 2, 3.$$

and λ3 σ∗ 2 � � ¼ �c4 σ∗ 2 � � 6¼ 0: We summarise the discussion above in the following theorem. Theorem 11. Suppose that the equilibrium point

$$E\_3(s^\*, \mathfrak{x}^\*, y^\*; \sigma) = E\_3\left(s^\*, \frac{1 - s^\*}{f\_1(s^\*)}, \frac{(2\mathfrak{x}^\* + 1)f\_2(\mathfrak{x}^\*) - \mathfrak{x}^\* - 1}{bD\_1(1 - f\_2(\mathfrak{x}^\*))}; \sigma\right)$$

exists and those assumptions similar as in (2) for the system (3) together with hypotheses H4 until H6 hold. Then the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating from E3 s∗; x∗; y∗ ð Þ ; σ for some values of σ in the neighbourhood of σ∗ 2:

#### 2.6. Discussion

We have proposed and analysed a simple nutrient-predator-prey model in a chemostat with intratrophic predation. This system consisted of the nutrient s, prey organisms x and predator organisms y. We conclude that intratrophic predation denoted as b does not affected the local stability and existence of E1 and E2 when no predator organisms involved. Next, we have shown that the washout equilibrium point E1ð Þ 1; 0; 0 (no prey and predator organisms present) is locally asymptotically stable if λ3 ¼ f 1ð Þ1 D1 � 1 < 0 holds. For stability of E2 and E3 of the system (3), some sufficient criteria or conditions are derived and satisfied. The points E2 and E3 are said to be asymptotically stable if all of their eigenvalues are less than zero.

[6] Pitchford J, Brindley J. Intratrophic predation in simple predator-prey models. Bulletin of

Stability and Hopf Bifurcation Analysis of a Simple Nutrient-Prey-Predator Model with Intratrophic Predation…

http://dx.doi.org/10.5772/intechopen.71624

67

[7] Mada Sanjaya WS, Salleh Z, Mamat M. Mathematical model of three species food chain with Holling Type-III functional response. International Journal of Pure and Applied

[8] Tee LS, Salleh Z. Hopf bifurcation of a nonlinear system derived from Lorenz system using normal form theory. International Journal of Applied Mathematics and Statistics.

[9] Hassard BD, Kazarinoff ND, Wan YH. Theory and Applications of Hopf Bifurcation.

[10] Freedman HI, Ruan S. Hopf bifurcation in three-species food chain models with group

[11] Rana SMS, Nasrin F, Hossain MI. Dynamics of a predator-prey interaction in chemostat with variable yield. Journal of Sustainability Science and Management. 2015;10(2):16-23

[12] Solomon ME. The natural control of animal populations. The Journal of Animal Ecology.

[14] Robinson CR. An Introduction to Dynamical Systems: Continuous and Discrete. Pearson

[15] Hirsch MW, Smale S, Devaney RL. Differential Equations, Dynamical Systems & an Introduction to Chaos. 2nd ed. Elsevier Academic Press, San Diego California: Elsevier; 2004

[16] Freedman HI, Mathsen RM. Persistence in predator-prey systems with ratio-dependent

[17] Butler GJ, Hsu SB, Waltman P. Coexistence of competing predators in a chemostat. Jour-

[18] Nani F, Freedman HI. A mathematical model of cancer treatment by immunotherapy.

[19] Freedman HI, Waltman P. Persistence in models of three interacting predator-prey

[20] Freedman HI, Ruan S. Uniform persistence in functional differential equations. Mathe-

[22] Tee LS, Salleh Z. Hopf bifurcation analysis of Zhou system. International Journal of Pure

[23] Marsden JE, McKracken M. The Hopf Bifurcation and its Applications. New York:

[21] Gantmacher FR. The Theory of Matrices. USA: Chelsea Publishing Company; 1964

predator influence. Bulletin of Mathematical Biology. 1993;55(4):817-827

[13] Hale JK, Kocak H. Dynamics and Bifurcations. New York: Springer; 1991

Prentice Hall, New Jersey: Pearson Education, Inc; 2004

nal of Mathematical Biology. 1983;17:133-151

Mathematical Biosciences. 2000;163:159-199

matical Biosciences. 1995;115:173-192

Springer-Verlag; 1976

and Applied Mathematics. 2015;104(1):1-18

populations. Mathematical Biosciences. 1984;68:213-231

Mathematical Biology. 1998;60:937-953; bu980053

USA: Cambridge University Press; 1981

defense. Mathematical Biosciences. 1992;111:73-87

2016;55(3):122-132

1949;18:1-32

Mathematics. 2013;89(5):647-657. DOI: 10.12737/ijpam.v89i5.1

In particular, we investigated the global stability analysis for E2: A suitable Lyapunov function L is defined and E2 is globally asymptotically stable if and only if the conditions in Theorem 8 holds. Global stability of E2 indicates that predator organisms y might be washout in the chemostat despites the initial prey and predator organisms' density levels. Next, in the study of the existence of positive interior equilibrium point E3, we presented some results regarding uniform persistence analysis. It has shown that the system (3) is uniformly persistence and thus, the positive interior equilibrium point E3 exists.

In the analysis of occurrence of Hopf bifurcation, Hopf bifurcation theorems in Hassard et al. [9] are applied. We have shown that the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating from E2 ζs; ð Þ 1�ζs D1 ; 0; σ and E3 s∗; x∗; y∗ ð Þ ; σ for some values of σ in the neighbourhoods of σ∗ 1 and σ∗ 2, respectively. The method used to obtain these results is similar to the method used in [3, 10].

#### Author details

Zabidin Salleh\* and Liyana Abd Rahim

\*Address all correspondence to<EMAIL_ADDRESS>
School of Informatics and Applied Mathematics, Universiti Malaysia Terengganu, Kuala Nerus, Terengganu, Malaysia

#### References

[6] Pitchford J, Brindley J. Intratrophic predation in simple predator-prey models. Bulletin of Mathematical Biology. 1998;60:937-953; bu980053

shown that the washout equilibrium point E1ð Þ 1; 0; 0 (no prey and predator organisms present)

system (3), some sufficient criteria or conditions are derived and satisfied. The points E2 and

In particular, we investigated the global stability analysis for E2: A suitable Lyapunov function L is defined and E2 is globally asymptotically stable if and only if the conditions in Theorem 8 holds. Global stability of E2 indicates that predator organisms y might be washout in the chemostat despites the initial prey and predator organisms' density levels. Next, in the study of the existence of positive interior equilibrium point E3, we presented some results regarding uniform persistence analysis. It has shown that the system (3) is uniformly persistence and

In the analysis of occurrence of Hopf bifurcation, Hopf bifurcation theorems in Hassard et al. [9] are applied. We have shown that the system (8) undergoes Hopf bifurcation in the first

E3 are said to be asymptotically stable if all of their eigenvalues are less than zero.

octant, which leads to a family of periodic solutions bifurcating from E2 ζs;

School of Informatics and Applied Mathematics, Universiti Malaysia Terengganu,

tion. New York: Cambridge University Press; 1995. pp. 1-7

[1] Smith HL, Waltman P. The Theory of the Chemostat: Dynamics of Microbial Competi-

[2] Li B, Kuang Y. Simple food chain in a chemostat with distinct removal rates. Journal of Mathematical Analysis and Applications. 2000;242:75-92. DOI: 10.1006 jmaa.1999.6655 [3] El-Sheikh MMA, Mahrouf SAA. Stability and bifurcation of a simple food chain in a chemostat with removal rates. Chaos, Solitons & Fractals. 2005;23:1475-1489. DOI: 10.1016/j.

[4] Hamra MA, Yadi K. Asymptotic behavior of a chemostat model with constant recycle sludge concentration. Acta Biotheoretica. 2017;65:233-252. DOI: 10.1007/s10441-017-9309-4

[5] Yang J, Tang G, Tang S. Modelling the regulatory system of a chemostat model with a threshold window. Mathematics and Computers in Simulation. 2017;132:220-235. DOI:

method used to obtain these results is similar to the method used in [3, 10].

E3 s∗; x∗; y∗ ð Þ ; σ for some values of σ in the neighbourhoods of σ∗

D1 � 1 < 0 holds. For stability of E2 and E3 of the

ð Þ 1�ζs D1 ; 0; σ

2, respectively. The

1 and σ∗

and

is locally asymptotically stable if λ3 ¼ f 1ð Þ1

66 Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

thus, the positive interior equilibrium point E3 exists.

Author details

References

Zabidin Salleh\* and Liyana Abd Rahim

Kuala Nerus, Terengganu, Malaysia

chaos.2004.06.079

10.1016/j.matcom.2016.08.005

\*Address all correspondence to<EMAIL_ADDRESS>

**Chapter 4**

Provisional chapter

**Sensitivity Analysis: A Useful Tool for Bifurcation**

DOI: 10.5772/intechopen.72345

Sensitivity analysis and bifurcation analysis are closely related to each other. In sensitivity analysis, especially global sensitivity analysis the effects of input parameter spaces on output quantities of interest are studied. On the other hand, in bifurcation analysis the critical points within feasible regions of parameters are detected where the long-term dynamics changes qualitatively. Prior to bifurcation analysis, it is important to identify the bifurcation parameters. In complex and computationally expensive problems which consist plenty of uncertain parameters, it is essential to find a set of bifurcation parameters before bifurcation analysis. Global sensitivity analysis is a powerful tool to identify the bifurcation parameters which contribute most on output uncertainty. Global sensitivity analysis is the first step toward bifurcation analysis which helps in dimension reduction during bifurcation analysis. As an example, in this chapter, a multi compartment, lumpedparameter model of an arm artery is considered and global sensitivity analysis (Sobol's

method) is applied to identify the bifurcation parameters of the arm arteries.

analysis, bifurcation parameters, Sobol's method

Keywords: lumped parameter model, arm arteries, sensitivity analysis, bifurcation

Sensitivity analysis and bifurcation analysis are closely related to each other. In sensitivity analysis, we study how the uncertainty in the output of a mathematical model or system (numerical or otherwise) can be apportioned to different sources of uncertainty in its inputs [1]. On the other hand, in bifurcation analysis the critical points within the feasible regions of parameters are detected where the long-term dynamics changes qualitatively [2]. Prior to the bifurcation analysis, it is important to identify the bifurcation parameters in complex and

> © 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

distribution, and reproduction in any medium, provided the original work is properly cited.

Sensitivity Analysis: A Useful Tool for Bifurcation

**Analysis**

Analysis

Raheem Gul and Stefan Bernhard

Raheem Gul and Stefan Bernhard

http://dx.doi.org/10.5772/intechopen.72345

Abstract

1. Introduction

Additional information is available at the end of the chapter

Additional information is available at the end of the chapter

#### **Sensitivity Analysis: A Useful Tool for Bifurcation Analysis** Sensitivity Analysis: A Useful Tool for Bifurcation Analysis

DOI: 10.5772/intechopen.72345

Raheem Gul and Stefan Bernhard Raheem Gul and Stefan Bernhard

Additional information is available at the end of the chapter Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/intechopen.72345

#### Abstract

Sensitivity analysis and bifurcation analysis are closely related to each other. In sensitivity analysis, especially global sensitivity analysis the effects of input parameter spaces on output quantities of interest are studied. On the other hand, in bifurcation analysis the critical points within feasible regions of parameters are detected where the long-term dynamics changes qualitatively. Prior to bifurcation analysis, it is important to identify the bifurcation parameters. In complex and computationally expensive problems which consist plenty of uncertain parameters, it is essential to find a set of bifurcation parameters before bifurcation analysis. Global sensitivity analysis is a powerful tool to identify the bifurcation parameters which contribute most on output uncertainty. Global sensitivity analysis is the first step toward bifurcation analysis which helps in dimension reduction during bifurcation analysis. As an example, in this chapter, a multi compartment, lumpedparameter model of an arm artery is considered and global sensitivity analysis (Sobol's method) is applied to identify the bifurcation parameters of the arm arteries.

Keywords: lumped parameter model, arm arteries, sensitivity analysis, bifurcation analysis, bifurcation parameters, Sobol's method

#### 1. Introduction

Sensitivity analysis and bifurcation analysis are closely related to each other. In sensitivity analysis, we study how the uncertainty in the output of a mathematical model or system (numerical or otherwise) can be apportioned to different sources of uncertainty in its inputs [1]. On the other hand, in bifurcation analysis the critical points within the feasible regions of parameters are detected where the long-term dynamics changes qualitatively [2]. Prior to the bifurcation analysis, it is important to identify the bifurcation parameters in complex and

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited. © 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

computationally expensive problems that consist plenty of uncertain parameters. Sensitivity analysis is a powerful tool to identify the bifurcation parameters which contribute most on output uncertainty. Also, sensitivity analysis helps in dimension reduction during the bifurcation analysis by fixing less influential parameters on their nominal values.

In this chapter, the main questions of interest are:

2. Lumped-parameter model of the arm arteries

Pressure and flow equations at non-terminal nodes:

q\_

q\_

p\_i ¼ qi � qiþ1 Ci

p\_ 6 ¼ q6 � q11 � q7 C6

Pressure and flow equations at terminal nodes:

q\_i ¼ pi�1 � pi � Riqi Li

11 ¼ p6 � p11 � R11q11 L11

13 ¼ p11 � p13 � R13q13 L13

Flow equation:

Pressure equation:

1. How to identify the bifurcation parameters in a model having plenty of input parameters? 2. Which parameters could be exempted from the bifurcation analysis (dimension reduction)? This chapter seeks to answers the above-mentioned questions using a simplified 5-steps procedure of uncertainty and sensitivity analysis. As an example, a multi-compartment, lumpedparameter model of arm arteries is considered [4] and global sensitivity analysis (Sobol's method) is applied to identify the bifurcation parameters (electrical) of the arm arteries.

In this section, the major arteries of the arm are divided in to number of non-terminal and terminal arterial segments (nodes). The total number of arterial segments, Ns ¼ 15 including 12 non-terminal and 3 terminal segments. Each non-terminal and terminal arterial segment is

Applying Kirchhoff's current and voltage laws on electrical representation of arm arteries, the

, i ¼ 1, 2, 3, …15 and i 6¼ 6, 11

, p\_ 11 ¼ q11 � q12 � q13

q\_in ¼ 2pin � 2pi � Riqin Li

q\_out ¼ 2pi � 2pout � 2Rbqout Li

p\_i ¼ qin � qout Ci

C11

, i ¼ 1, 2, 3, …15 and i 6¼ 11, 13

, i ¼ 10, 12, 15

ð Þ at bifurcation

Sensitivity Analysis: A Useful Tool for Bifurcation Analysis

http://dx.doi.org/10.5772/intechopen.72345

71

(1)

(2)

(3)

represented by its corresponding non-terminal and terminal electrical circuit.

following mathematical equations for pressure and flow are obtained:

Sensitivity analysis can be divided into two categories, local sensitivity analysis (LSA) and global sensitivity analysis (GSA). In LSA a parameter value is perturbed around its nominal values at a time, keeping other parameters fixed on their nominal values. The procedure is repeated for all parameters one by one to study their impact on output variables. LSA techniques are simple, easy to implement and computationally less expensive. On the other hand, LSA is not suitable for non-linear models and does not explore the impact of entire parameter spaces and their interactions effects on output variables [3, 4]. In order to overcome the limitations associated with the LSA, GSA can be used. In GSA, the analysis is performed over entire feasible regions of the input parameters and quantifies the impact of parameter interactions on output variables. The only deficiency related to the GSA is its computational cost [5–12] (Figure 1).

Figure 1. A simplified 5-step procedure to identify the bifurcation parameters using global sensitivity analysis.

In this chapter, the main questions of interest are:

1. How to identify the bifurcation parameters in a model having plenty of input parameters?

2. Which parameters could be exempted from the bifurcation analysis (dimension reduction)?

This chapter seeks to answers the above-mentioned questions using a simplified 5-steps procedure of uncertainty and sensitivity analysis. As an example, a multi-compartment, lumpedparameter model of arm arteries is considered [4] and global sensitivity analysis (Sobol's method) is applied to identify the bifurcation parameters (electrical) of the arm arteries.

#### 2. Lumped-parameter model of the arm arteries

In this section, the major arteries of the arm are divided in to number of non-terminal and terminal arterial segments (nodes). The total number of arterial segments, Ns ¼ 15 including 12 non-terminal and 3 terminal segments. Each non-terminal and terminal arterial segment is represented by its corresponding non-terminal and terminal electrical circuit.

Applying Kirchhoff's current and voltage laws on electrical representation of arm arteries, the following mathematical equations for pressure and flow are obtained:

Pressure and flow equations at non-terminal nodes:

Flow equation:

computationally expensive problems that consist plenty of uncertain parameters. Sensitivity analysis is a powerful tool to identify the bifurcation parameters which contribute most on output uncertainty. Also, sensitivity analysis helps in dimension reduction during the bifurca-

Sensitivity analysis can be divided into two categories, local sensitivity analysis (LSA) and global sensitivity analysis (GSA). In LSA a parameter value is perturbed around its nominal values at a time, keeping other parameters fixed on their nominal values. The procedure is repeated for all parameters one by one to study their impact on output variables. LSA techniques are simple, easy to implement and computationally less expensive. On the other hand, LSA is not suitable for non-linear models and does not explore the impact of entire parameter spaces and their interactions effects on output variables [3, 4]. In order to overcome the limitations associated with the LSA, GSA can be used. In GSA, the analysis is performed over entire feasible regions of the input parameters and quantifies the impact of parameter interactions on output variables. The only deficiency related to the GSA is its computational cost [5–12]

Figure 1. A simplified 5-step procedure to identify the bifurcation parameters using global sensitivity analysis.

tion analysis by fixing less influential parameters on their nominal values.

70 Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

(Figure 1).

$$\begin{aligned} \dot{q}\_{i} &= \frac{p\_{i-1} - p\_i - R\_i q\_i}{L\_i}, \quad i = 1, 2, 3, \dots 15 \text{ and } i \neq 11, 13\\ \dot{q}\_{11} &= \frac{p\_6 - p\_{11} - R\_{11} q\_{11}}{L\_{11}}\\ \dot{q}\_{13} &= \frac{p\_{11} - p\_{13} - R\_{13} q\_{13}}{L\_{13}} \end{aligned} \tag{1}$$

Pressure equation:

$$\begin{aligned} \dot{p}\_i &= \frac{q\_i - q\_{i+1}}{\mathcal{C}\_i}, \quad i = 1, 2, 3, \dots 15 \text{ and } i \neq 6, 11\\ \dot{p}\_6 &= \frac{q\_6 - q\_{11} - q\_7}{\mathcal{C}\_6}, \quad \dot{p}\_{11} = \frac{q\_{11} - q\_{12} - q\_{13}}{\mathcal{C}\_{11}} \end{aligned} \tag{2}$$

Pressure and flow equations at terminal nodes:

$$\begin{aligned} \dot{q}\_{in} &= \frac{2p\_{in} - 2p\_i - R\_i q\_{in}}{L\_i} \\ \dot{p}\_i &= \frac{q\_{in} - q\_{out}}{C\_i} \\ \dot{q}\_{out} &= \frac{2p\_i - 2p\_{out} - 2R\_b q\_{out}}{L\_i}, \quad i = 10, 12, 15 \end{aligned} \tag{3}$$

The value of boundary resistance ð Þ Rb on three terminal nodes is 3:24 � 109 kgs�1 m�4, ν ¼ 0:004 kgs�1 m�1 and ρ ¼ 1050 kgm�3 [4, 5, 13, 14].

Table 1. Numerical values of parameters for each node of the arm arteries (shown in Figure 2).

where, Ri, Ci and Li is the blood flow resistance, compliance of the vessel and blood inertia of i th segment of the arm arteries respectively. The electrical parameters Ri; Ci ð Þ ; Li of i th segments are related with structural parameters Ei; li; di ð Þ ; hi as,

$$R\_i = \frac{8\nu l\_i}{\pi \left(\frac{d}{2}\right)^{4'}}, \quad \mathcal{C}\_i = \frac{\rho l\_i}{\pi \left(\frac{d}{2}\right)^2}, \quad L\_i = \frac{2\pi \left(\frac{d}{2}\right)^2 l\_i}{E\_i l\_i} \tag{4}$$

outputs caused by uncertainty of its inputs. Whereas, sensitivity analysis study the impact of input quantities of interest (QoI) on output quantities of interest (QoI). In this study, the input (QoI) are electrical parameters (Ri, Ci, Li) and output (QoI) are pressure and flow at each node of the arm arteries. Further, for uncertainty analysis Latin hypercube sampling (LHS) is used and variance-decomposition method (Sobol's method) is used for global sensitivity analysis (GSA). Compared to the high-dimensional cardiovascular models (3D, 2D, 1D), lumped-parameter models of the cardiovascular system (CVS) are computationally less expensive, therefore they are suitable for GSA. In our previous studies, we found that for lumped-parameter models of

Figure 2. Model geometry of arm artery (A), with total number of arterial segments, Ns = 15, including 12 non-terminal and 3 terminal segments. Each non-terminal and terminal segment is represented by its corresponding non-terminal (C) and terminal electrical circuits (D). Pressure waves is used as an input boundary condition (B) and pout = 15 mmHg which is mean venous pressure used to calculate boundary outflow. The parameter values of each arterial segment are given in Table 1.

Sensitivity Analysis: A Useful Tool for Bifurcation Analysis

http://dx.doi.org/10.5772/intechopen.72345

73

where, Ei is the Young modulus, li denotes length of the vessel, d is the diameter of the vessel and hi represents the wall thickness of i th segment of the vessel. Moreover, ν (0.004 Pa s) is the blood viscosity and ρ (1050 kgm�3 ) is the blood density. The nominal values of all parameters of arm segments are given in Table 1. The geometry along with the values of the parameters is taken from [13, 14].

#### 3. Uncertainty and sensitivity analysis

Uncertainty analysis (UA) and sensitivity analysis (SA) are closely related; however they represent two different disciplines. Uncertainty analysis assesses the uncertainty in model

Figure 2. Model geometry of arm artery (A), with total number of arterial segments, Ns = 15, including 12 non-terminal and 3 terminal segments. Each non-terminal and terminal segment is represented by its corresponding non-terminal (C) and terminal electrical circuits (D). Pressure waves is used as an input boundary condition (B) and pout = 15 mmHg which is mean venous pressure used to calculate boundary outflow. The parameter values of each arterial segment are given in Table 1.

where, Ri, Ci and Li is the blood flow resistance, compliance of the vessel and blood inertia of

π d 2

where, Ei is the Young modulus, li denotes length of the vessel, d is the diameter of the vessel

of arm segments are given in Table 1. The geometry along with the values of the parameters is

Uncertainty analysis (UA) and sensitivity analysis (SA) are closely related; however they represent two different disciplines. Uncertainty analysis assesses the uncertainty in model

2 , Li ¼ 2π d

2 2 li Eili

m�4 � 106 kg�1

s2m4 � 10�11 kgm�4 � 106

th segment of the vessel. Moreover, ν (0.004 Pa s) is the

) is the blood density. The nominal values of all parameters

th segments

m�1 and

m�4, ν ¼ 0:004 kgs�1

(4)

th segment of the arm arteries respectively. The electrical parameters Ri; Ci ð Þ ; Li of i

Nodes E l d h R C L

1 4 6.1 7.28 6.2 3.539 7.454 1.539 2 4 5.6 6.28 5.7 5.868 4.778 1.898 3 4 6.3 5.64 5.5 10.15 4.035 2.648 4 4 6.3 5.32 5.3 12.82 3.514 2.976 5 4 6.3 5 5.2 16.43 2.974 3.369 6 4 4.6 4.72 5 15.10 1.9 2.76 7 8 7.1 3.48 4.4 78.90 0.667 7.838 8 8 7.1 3.24 4.3 105 0.531 9.042 9 8 7.1 3 4.2 142.9 0.448 10.55 10 8 2 2.84 4.1 55.11 0.1207 3.647 11 8 2 4.3 4.9 31.94 1.067 4.844 12 16 6.7 1.82 2.8 1173 0.0834 31.88 13 8 7.9 4.06 4.7 40.19 0.9366 5.434 14 8 6.7 3.48 4.6 50.22 0.80 6.075 15 8 3.7 3.66 4.5 33.60 0.3958 3.693

s�2 �105 m �10�2 m �10�3 m �10�4 kgs�1

72 Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

4 , Ci ¼ ρli

are related with structural parameters Ei; li; di ð Þ ; hi as,

and hi represents the wall thickness of i

3. Uncertainty and sensitivity analysis

blood viscosity and ρ (1050 kgm�3

taken from [13, 14].

ρ ¼ 1050 kgm�3 [4, 5, 13, 14].

Ri ¼ 8νli π d 2

The value of boundary resistance ð Þ Rb on three terminal nodes is 3:24 � 109 kgs�1

Table 1. Numerical values of parameters for each node of the arm arteries (shown in Figure 2).

i

units kgm�2

outputs caused by uncertainty of its inputs. Whereas, sensitivity analysis study the impact of input quantities of interest (QoI) on output quantities of interest (QoI). In this study, the input (QoI) are electrical parameters (Ri, Ci, Li) and output (QoI) are pressure and flow at each node of the arm arteries. Further, for uncertainty analysis Latin hypercube sampling (LHS) is used and variance-decomposition method (Sobol's method) is used for global sensitivity analysis (GSA).

Compared to the high-dimensional cardiovascular models (3D, 2D, 1D), lumped-parameter models of the cardiovascular system (CVS) are computationally less expensive, therefore they are suitable for GSA. In our previous studies, we found that for lumped-parameter models of the CVS, the Sobol's method is computationally less expensive as compared to the other variance-decomposition methods, like FAST and sparse grid stochastic collocation method based on Smolyak algorithm [5, 12].

where, V is the variance operator. The relationship between functions and partial variance are

ð Þ xi � �

� � � � � Vi � V ¼ V f ij xi; xj

� � � �

Sensitivity Analysis: A Useful Tool for Bifurcation Analysis

http://dx.doi.org/10.5772/intechopen.72345

� � þ … þ S1,2,3,…,Kð Þ x1; x2; x3;…; xK (10)

ð Þ Yjx � i � �

(9)

75

(11)

th parameter on output

th parameters on output uncertainty.

V (12)

th parameter means main effect

th parameter. In this study, the interaction

ð Þ¼ Ex�ið Þ Yjxi V f i

Vij ¼ Vxi,xj Ex�ij Yjxi; xj

X k

X k

Sij xi; xj

Si ¼ Vij

Sij ¼ Vij V

V , and

th and j

V ¼ 1 � Vx�i Exi

In general, the main effect is used identify the most influential parameters (bifurcation parameters) and the total effect is taken into account for those parameters which are exempted from

effects of parameters on the output (QoI) are negligible, therefore the main effects are used for

In this section, a detailed working algorithm is presented to compute the main effect, Si using

1. Generate a random numbers matrix of row dimension 2K and column length N (the sample size) and split into two independent sampling matrices, Að Þ N;K and Bð Þ N;K by

i<j

where, Si is the main effect (first order sensitivity index) of the i

ð Þ Yjx � i � �

i

given by,

Where,

Vi ¼ Vxi

…

Dividing both sides of the Eq. (8) by V, we get:

Sið Þþ xi

uncertainty and Sij is the interaction effect of i

Further, the total sensitivity index, STi can be calculated as,

STi ¼ Ex�i Vxi

bifurcation analysis (factor fixing). The total effect, STi of the i

the Monte Carlo simulations, we follow the steps, given in [1, 15].

using LHS. Where, K is the number of uncertain model parameters.

plus higher-order effect due to interactions of the i

factor fixing and ranking of bifurcation parameters.

3.2. Algorithm to compute sensitivity indices

1 ¼ X k

i

#### 3.1. The method of Sobol

The method of Sobol is the variance-decomposition method used for global sensitivity analysis. The method decomposes the output variance of a system or model into fractions and assigns them to the inputs factors. For example, given a model of the form Y ¼ f Xð Þ¼ f xð Þ 1; x2;…; xk , where X is the vector of K uncertain parameters, which are independently generated within a unit hypercube i.e. xi ∈½ � 0; 1 k for i ¼ 1, 2, 3,…, K. Compared to the other GSA methods, the Sobol's method is one of the most commonly used variance-decomposition method, because of its ease of implementation. The method is primarily based on the decomposition of output Y into summands of elementary functions in terms of increasing dimensionality [1, 8],

$$f(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_k) = f\_0 + \sum\_{i}^k f\_i(\mathbf{x}\_i) + \sum\_{i}^k \sum\_{i$$

In Eq. (5), f is integrable, f 0 is a constant, f i is a function of xi, f ij is a function of xi and xj and so on. Furthermore, all the terms in the functional decomposition are orthogonal, which leads toward the following definitions of the terms of the functional decomposition in term of conditional expected values.

$$\begin{aligned} f\_0 &= E(Y) \\ f\_i(\mathbf{x}\_i) &= E\_{\mathbf{x} \sim i}(Y|\mathbf{x}\_i) - f\_0 \\ f\_{ij}(\mathbf{x}\_i, \mathbf{x}\_j) &= E\_{\mathbf{x} \sim \vec{\eta}}(Y|\mathbf{x}\_i, \mathbf{x}\_j) - f\_0 - f\_i - f\_j \\ \dots & \dots \end{aligned} \tag{6}$$

where, E describes the mathematical expectation and x�i denotes all parameters except xi and so on. The total unconditional variance can be obtained by,

$$V = \int\_{\Omega} f^2(X)dx - f\_0^2 \tag{7}$$

From Eq. (7), the total unconditional variance can be decomposed in a similar manner like in Eq. (5) as,

$$V = \sum\_{i}^{k} V\_{i}(\mathbf{x}\_{i}) + \sum\_{i}^{k} \sum\_{i$$

where, V is the variance operator. The relationship between functions and partial variance are given by,

$$\begin{aligned} V\_i &= V\_{\mathbf{x}\_i}(\mathbb{E}\_{\mathbf{x}\sim i}(Y|\mathbf{x}\_i)) = V\left(f\_i(\mathbf{x}\_i)\right) \\ V\_{\vec{\eta}} &= V\_{\mathbf{x}\_i\mathbf{x}\_{\vec{\eta}}}\left(\mathbb{E}\_{\mathbf{x}\sim \vec{\eta}}\left(Y|\mathbf{x}\_i, \mathbf{x}\_{\vec{\eta}}\right)\right) - V\_i - V = V\left(f\_{\vec{\eta}}(\mathbf{x}\_i, \mathbf{x}\_{\vec{\eta}})\right) \\ &\dots \end{aligned} \tag{9}$$

Dividing both sides of the Eq. (8) by V, we get:

$$\mathbf{1} = \sum\_{i}^{k} \mathbf{S}\_{i}(\mathbf{x}\_{i}) + \sum\_{i}^{k} \sum\_{i$$

Where,

the CVS, the Sobol's method is computationally less expensive as compared to the other variance-decomposition methods, like FAST and sparse grid stochastic collocation method

The method of Sobol is the variance-decomposition method used for global sensitivity analysis. The method decomposes the output variance of a system or model into fractions and assigns them to the inputs factors. For example, given a model of the form Y ¼ f Xð Þ¼ f xð Þ 1; x2;…; xk , where X is the vector of K uncertain parameters, which are independently generated within a unit hypercube i.e. xi ∈½ � 0; 1 k for i ¼ 1, 2, 3,…, K. Compared to the other GSA methods, the Sobol's method is one of the most commonly used variance-decomposition method, because of its ease of implementation.  Furthermore, all the terms in the functional decomposition are orthogonal, which leads toward the following definitions of the terms of the functional decomposition in term of

where, E describes the mathematical expectation and x�i denotes all parameters except xi and

From Eq. (7), the total unconditional variance can be decomposed in a similar manner like in

ð Þ X dx � f

i

ð Þ¼ xi Ex�ið Þ� Yjxi f 0

V ¼ ð

ΩK f 2

Vij xi; xj

� � ¼ Ex�ij Yjxi; xj

i f i ð Þþ xi

f 0 ¼ E Yð Þ

f ij xi; xj

f i

…

so on. The total unconditional variance can be obtained by,

3.1. The method of Sobol

ality [1, 8],

Eq. (5) as,

$$\begin{aligned} S\_i &= \frac{V\_{\vec{ij}}}{V}, \text{ and} \\ S\_{\vec{ij}} &= \frac{V\_{\vec{ij}}}{V} \end{aligned} \tag{11}$$

where, Si is the main effect (first order sensitivity index) of the i th parameter on output uncertainty and Sij is the interaction effect of i th and j th parameters on output uncertainty. Further, the total sensitivity index, STi can be calculated as,

$$\mathbf{S}\_{T\_i} = \frac{E\_{\mathbf{x} \sim i} \left( V\_{\mathbf{x}\_i}(\mathbf{Y}|\mathbf{x} \sim i) \right)}{V} = 1 - \frac{V\_{\mathbf{x} \sim i} \left( E\_{\mathbf{x}\_i}(\mathbf{Y}|\mathbf{x} \sim i) \right)}{V} \tag{12}$$

In general, the main effect is used identify the most influential parameters (bifurcation parameters) and the total effect is taken into account for those parameters which are exempted from bifurcation analysis (factor fixing). The total effect, STi of the i th parameter means main effect plus higher-order effect due to interactions of the i th parameter. In this study, the interaction effects of parameters on the output (QoI) are negligible, therefore the main effects are used for factor fixing and ranking of bifurcation parameters.

#### 3.2. Algorithm to compute sensitivity indices

In this section, a detailed working algorithm is presented to compute the main effect, Si using the Monte Carlo simulations, we follow the steps, given in [1, 15].

1. Generate a random numbers matrix of row dimension 2K and column length N (the sample size) and split into two independent sampling matrices, Að Þ N;K and Bð Þ N;K by using LHS. Where, K is the number of uncertain model parameters.

$$A(\mathbf{N}, K) = \begin{bmatrix} \mathbf{x}\_{11} & \mathbf{x}\_{12} & \dots & \mathbf{x}\_{1K} \\ \mathbf{x}\_{21} & \mathbf{x}\_{22} & \dots & \mathbf{x}\_{2K} \\ \dots & \dots & \dots & \dots \\ \mathbf{x}\_{N1} & \mathbf{x}\_{N2} & \dots & \mathbf{x}\_{NN} \end{bmatrix} \tag{13}$$

$$\begin{bmatrix} \mathbf{x}\_{1(K+1)} & \mathbf{x}\_{1(K+2)} & \dots & \mathbf{x}\_{1(2K)} \end{bmatrix}$$

5. Finally, the main effect, Si of each parameter on the state variables is calculated.

X Nt

ð Þ t; j; t , i ¼ 1, 2, …, K (19)

Sensitivity Analysis: A Useful Tool for Bifurcation Analysis

http://dx.doi.org/10.5772/intechopen.72345

77

t¼0 Sti

In Eq. (19), NTs is the number of output variables (pressure and flow time series at all

The results of the UA and SA are greatly affected by the choice of input parameters distributions. In principle, the parameters distributions should be estimated using medical data. Unfortunately, the medical data is not easy to obtained. The input parameters distributions could be chosen according to the expert opinion or using the data from the literature. Due to limited data availability, here in this work the input parameters are randomized within �10%

The method of Sobol requires Nð Þ K þ 2 number of model simulations to compute Si. The main effect, Si is computed for N = [500, 1000, 2000, 3000, 4000] model runs. It is observed that, when the total number of simulations run N increases from 3000 then the sensitivity indices (Si) become stable [18]. Therefore, the minimum number of simulations for each parameter to

In this section, the sensitivity results based on main effect Si are presented. In order to calculate sensitivity time series, the method of Sobol is applied on each time point of the output QoI i.e. pressure and flow waves at each location of the arm arteries. For each parameter, there are two sensitivity time series at each segment of the arm arteries, one for the pressure and one for the flow. In total, K � NTs ¼ 45 � 33 ¼ 1485 sensitivity time series are obtained. In order to represent the sensitivity results in a compact way, mean absolute values of each pressure and flow sensitivity time series per parameter is taken. In this way, a matrix of dimension 45 � 33 is acquired, where each entry of the matrix represents the mean absolute values of pressure and flow sensitivity time series per parameter, see Figure 3. The numbers in the boxes show the impact (%) on the output (pressure and flow) when input parameters Ri;Ci ð Þ ; Li are randomized within the feasible ranges of �10%. The parameters having main effect, Si > 10% on output QoI are not shown in the Figure 3. Each row in Figure 3 represents the ranking of influential (bifurcation) parameters. For convenience, the electrical parameters Ri;Ci ð Þ ; Li , i ¼ 1, 2, 3, …, 15 that have impact greater than 10% on pressure and flow are considered as bifurcation parameters which further can be used in bifurcation analysis. For example, for pressure

j¼1

range of their base (nominal) values using Latin hypercube sampling (LHS).

1 Nt X NTs

Si ¼ 1 NTs

locations) and Nt is the number of time-points [12].

achieve convergence of sensitivity indices is around 3000.

3.3. Input parameters distribution

3.4. Convergence of sensitivity indices

4. Results and discussion

$$B(\mathbf{N}, K) = \begin{bmatrix} \mathbf{x}\_{1(\mathsf{K}+1)} & \mathbf{x}\_{1(\mathsf{K}+2)} & \dots & \mathbf{x}\_{1(\mathsf{K})} \\ \mathbf{x}\_{2(\mathsf{K}+1)} & \mathbf{x}\_{2(\mathsf{K}+2)} & \dots & \mathbf{x}\_{2(2\mathsf{K})} \\ \dots & \dots & \dots & \dots \\ \mathbf{x}\_{\mathsf{N}(\mathsf{K}+1)} & \mathbf{x}\_{\mathsf{N}(\mathsf{K}+2)} & \dots & \mathbf{x}\_{\mathsf{N}(2\mathsf{K})} \end{bmatrix} \tag{14}$$

2. Define matrix Ci, which is matrix A except the i th column of matrix B.

$$\mathbf{C}\_{i}(\mathbf{N},K) = \begin{bmatrix} \mathbf{x}\_{11} & \mathbf{x}\_{12} & \dots & \mathbf{x}\_{1}(K+i) & \dots & \mathbf{x}\_{1K} \\ \mathbf{x}\_{21} & \mathbf{x}\_{22} & \dots & \mathbf{x}\_{2}(K+i) & \dots & \mathbf{x}\_{2K} \\ \dots & \dots & \dots & \dots & \dots & \dots \\ \mathbf{x}\_{N1} & \mathbf{x}\_{N2} & \dots & \mathbf{x}\_{N}(K+i) & \dots & \mathbf{x}\_{NK} \end{bmatrix} \tag{15}$$

$$\begin{split} S\_{i} = \frac{V\_{i}}{V} = \frac{V\_{\text{xi}i}(E\_{\text{x}\sim i}(\mathbf{Y}|\mathbf{x}\_{i}))}{V} = \frac{V - \frac{1}{2\mathbf{N}}\sum\_{n=1}^{\text{N}} \left( Y\_{\text{B}}^{(n)} - Y\_{\text{C}\_{i}}^{(n)} \right)^{2}}{V} \\ = 1 - \frac{1}{2\mathbf{N} \times V} \sum\_{n=1}^{\text{N}} \left( Y\_{\text{B}}^{(n)} - Y\_{\text{C}\_{i}}^{(n)} \right)^{2} \end{split} \tag{16}$$

where,

$$V = \frac{1}{N} \sum\_{n=1}^{N} \left( Y\_{\mathcal{B}}^{(n)} \right)^{2} - E^{2} \tag{17}$$

and

$$E = \left(\frac{1}{\mathcal{N}} \sum\_{n=1}^{N} Y\_{\mathcal{B}}^{(n)}\right)^2 \tag{18}$$

The total variance ð Þ V and the expectation ð Þ E are also calculated at each time-point of pressure and flow waves with respect to each parameter.

5. Finally, the main effect, Si of each parameter on the state variables is calculated.

$$S\_{l} = \frac{1}{\mathbf{N}\_{T\_{s}}} \frac{1}{\mathbf{N}\_{l}} \sum\_{j=1}^{\mathbf{N}\_{l}} \sum\_{t=0}^{\mathbf{N}\_{l}} S\_{l\natural}(t, j, t), i = 1, 2, \dots, K \tag{19}$$

In Eq. (19), NTs is the number of output variables (pressure and flow time series at all locations) and Nt is the number of time-points [12].

#### 3.3. Input parameters distribution

Að Þ¼ N;K

Bð Þ¼ N;K

76 Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

V ¼ Vxið Þ Ex�ið Þ Yjxi

pressure and flow waves with respect to each parameter.

 Define matrix Ci, which is matrix A except the i

Cið Þ¼ N;K

model runs Nð Þ ¼ 4000 .

where,

and

estimator offered by Jansen [15–17].

 Compute and save model runs for all parameter spaces using matrices A, B and Ci i.e. YA t; Ts ð Þ¼ ; N f Að Þ, YB t; Ts ð Þ¼ ; N f Bð Þ and YCi t; Ts ð Þ¼ ; N;K f Cð Þi , where, t are the time points for one heart beat with period tp ¼ 0:8s, Ts represents the state variables (pressure and flow time series at six locations of arm artery Nð Þ Ts ¼ 15 and Nis the total number of

4. For the time dependent model outputs, we compute the time dependent main sensitivity index, of each parameter at each time-point of the pressure and flow waves, using the

> V � 1 2N P N n¼1

¼ 1 � 1

Yð Þ n B � �2

n¼1

The total variance ð Þ V and the expectation ð Þ E are also calculated at each time-point of

Yð Þ n B !2

2N � V

Yð Þ n B � Yð Þ n Ci � �2

> X N

Yð Þ n B � Yð Þ n Ci � �2

� E2 (17)

n¼1

V

th column of matrix B.

(13)

(14)

(15)

(16)

(18)

The results of the UA and SA are greatly affected by the choice of input parameters distributions. In principle, the parameters distributions should be estimated using medical data. Unfortunately, the medical data is not easy to obtained. The input parameters distributions could be chosen according to the expert opinion or using the data from the literature. Due to limited data availability, here in this work the input parameters are randomized within �10% range of their base (nominal) values using Latin hypercube sampling (LHS).

#### 3.4. Convergence of sensitivity indices

The method of Sobol requires Nð Þ K þ 2 number of model simulations to compute Si. The main effect, Si is computed for N = [500, 1000, 2000, 3000, 4000] model runs. It is observed that, when the total number of simulations run N increases from 3000 then the sensitivity indices (Si) become stable [18]. Therefore, the minimum number of simulations for each parameter to achieve convergence of sensitivity indices is around 3000.

#### 4. Results and discussion

In this section, the sensitivity results based on main effect Si are presented. In order to calculate sensitivity time series, the method of Sobol is applied on each time point of the output QoI i.e. pressure and flow waves at each location of the arm arteries. For each parameter, there are two sensitivity time series at each segment of the arm arteries, one for the pressure and one for the flow. In total, K � NTs ¼ 45 � 33 ¼ 1485 sensitivity time series are obtained. In order to represent the sensitivity results in a compact way, mean absolute values of each pressure and flow sensitivity time series per parameter is taken. In this way, a matrix of dimension 45 � 33 is acquired, where each entry of the matrix represents the mean absolute values of pressure and flow sensitivity time series per parameter, see Figure 3. The numbers in the boxes show the impact (%) on the output (pressure and flow) when input parameters Ri;Ci ð Þ ; Li are randomized within the feasible ranges of �10%. The parameters having main effect, Si > 10% on output QoI are not shown in the Figure 3. Each row in Figure 3 represents the ranking of influential (bifurcation) parameters. For convenience, the electrical parameters Ri;Ci ð Þ ; Li , i ¼ 1, 2, 3, …, 15 that have impact greater than 10% on pressure and flow are considered as bifurcation parameters which further can be used in bifurcation analysis. For example, for pressure

Figure 3. Main effect sensitivity of 45-electrical parameters Ri; Ci ð Þ ; Li on pressure and flow time series pi ; qi at 15-segments ð Þ N1; N2; …; N15 of the arm arteries with number of simulations run per parameter, N is 4000. In total N(K + 2) = 4000(45 + 2) = 188,000 = 0.188 million of simulations run are required to compute the main and total sensitivity indices. The total time taken to compute the sensitivity indices is approximately 3 hours.

at node-2, L1 and L2 are the bifurcation parameters, see in Figure 4 (top). Whereas, for flow at node -2, L1 and L2 are considered as bifurcation parameters, see Figure 4 (bottom).

In a similar fashion, each row of Figure 3 represents the ranking of bifurcation parameters which further can be used in bifurcation analysis. The parameters which have main effect Si < 10% can be exempted from the bifurcation analysis. The criteria for factor fixing vary from problem to problem.

#### 5. Conclusion

In this chapter, a 5-step procedure of global sensitivity analysis is presented to identify the bifurcation parameters in a lumped-parameter model of the arm arteries. Moreover, the proposed procedure can be applied on any morphology or structure of the systemic circulation (carotid bifurcation, aorta or complete systemic circulation). The results of sensitivity analysis are useful to identify and rank the bifurcation parameters, as well as help which parameters could be exempted from the bifurcation analysis. In this particular example of the arm arteries, 23 out of 45 parameters can be excluded from the bifurcation analysis. Whereas, 22 identified as bifurcation parameters, which further can be used/studied in the bifurcation analysis.

Figure 4. Ranking of bifurcation parameters Ri;Ci ð Þ ; Li in complete arm arteries for pressure (top) and flow (bottom) at node-2. It can be clearly seen that, L1, L2 and L1, L3 are considered as bifurcation parameters for pressure and flow at

Sensitivity Analysis: A Useful Tool for Bifurcation Analysis

http://dx.doi.org/10.5772/intechopen.72345

79

node-2 respectively.

at node-2, L1 and L2 are the bifurcation parameters, see in Figure 4 (top). Whereas, for flow at

15-segments ð Þ N1; N2; …; N15 of the arm arteries with number of simulations run per parameter, N is 4000. In total N(K + 2) = 4000(45 + 2) = 188,000 = 0.188 million of simulations run are required to compute the main and total sensitivity

Figure 3. Main effect sensitivity of 45-electrical parameters Ri; Ci ð Þ ; Li on pressure and flow time series pi

; qi at

In a similar fashion, each row of Figure 3 represents the ranking of bifurcation parameters which further can be used in bifurcation analysis. The parameters which have main effect Si < 10% can be exempted from the bifurcation analysis. The criteria for factor fixing vary

In this chapter, a 5-step procedure of global sensitivity analysis is presented to identify the bifurcation parameters in a lumped-parameter model of the arm arteries. Moreover, the proposed procedure can be applied on any morphology or structure of the systemic circulation (carotid bifurcation, aorta or complete systemic circulation). The results of sensitivity analysis are useful to identify and rank the bifurcation parameters, as well as help which parameters could be exempted from the bifurcation analysis. In this particular example of the arm arteries, 23 out of 45 parameters can be excluded from the bifurcation analysis. Whereas, 22 identified as bifurcation parameters, which further can be used/studied in the bifurcation analysis.

node -2, L1 and L2 are considered as bifurcation parameters, see Figure 4 (bottom).

indices. The total time taken to compute the sensitivity indices is approximately 3 hours.

78 Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

from problem to problem.

5. Conclusion

Figure 4. Ranking of bifurcation parameters Ri;Ci ð Þ ; Li in complete arm arteries for pressure (top) and flow (bottom) at node-2. It can be clearly seen that, L1, L2 and L1, L3 are considered as bifurcation parameters for pressure and flow at node-2 respectively.

#### Author details

Raheem Gul1 \* and Stefan Bernhard2

\*Address all correspondence to<EMAIL_ADDRESS>
1 COMSATS Institute of Information Technology, Abbottabad, Pakistan

2 Department of Electrical Engineering and Information Technology, Pforzheim University of Applied Sciences, Germany

[12] Gul R, Bernhard S. Parametric uncertainty and global sensitivity analysis in a model of the carotid bifurcation: Identification and ranking of most sensitive model parameters.

Sensitivity Analysis: A Useful Tool for Bifurcation Analysis

http://dx.doi.org/10.5772/intechopen.72345

81

[13] Noordergraaf A, Verdouw PD, Boom HBK. The use of an analog computer in a circula-

[14] Westerhof N, Noordergraaf A. Arterial viscoelasticity: A generalized model, effect of input impedance and wave travel in the systemic tree. Journal of Biomechanics. 1970;3:

[15] Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto F, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index.

[16] Jansen MJW. Analysis of variance designs for model output. Computer Physics Commu-

[17] Jansen MJW, Rossing WAH, Daamen RA. Monte Carlo estimation of uncertainty contributions from several independent multivariate sources. In: Predictability and Nonlinear Modeling in Natural Sciences and Economics. Dordrecht: Kluwer Academic Publishers;

[18] Yang J. Convergence and uncertainty analyses in Monte-Carlo based analysis. Environ-

tion model. Progress in Cardiovascular Diseases. 1963;5(5):419-439

Computer Physics Communications. 2010;181:259-270

mental Modeling and Software. 2011;26:444-457

Mathematical Biosciences. 2015;269:104-116

357-379

nications. 1999;117:35-43

1994. pp. 334-343

#### References

[12] Gul R, Bernhard S. Parametric uncertainty and global sensitivity analysis in a model of the carotid bifurcation: Identification and ranking of most sensitive model parameters. Mathematical Biosciences. 2015;269:104-116

Author details

Applied Sciences, Germany

& Sons Ltd; 2008

Biology. 2011;5(6):336-346

metrics. 1991;33:161-174

medical Engineering. 2013:1-24

\* and Stefan Bernhard2

\*Address all correspondence to<EMAIL_ADDRESS>
80 Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

Vol. 112. New York: Springer Verlag; 2004

1 COMSATS Institute of Information Technology, Abbottabad, Pakistan

2 Department of Electrical Engineering and Information Technology, Pforzheim University of

[1] Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, Saisana M, Tarantola S. Global Sensitivity Analysis. The Atrium, Chichester, West Sussex, England: John Wiley

[2] Kuznetsov YA. Elements of Applied Bifurcation Theory. Applied Mathematical Sciences.

[3] Zi Z. Sensitivity analysis approaches applied to systems biology models. IET Systems

[4] Gul R, Schütte, Bernhard S. Mathematical modeling and sensitivity analysis of arterial anastomosis in arm arteries. Applied Mathematical Modeling. 2016;40:7724-7738

[5] Gul R. Mathematical modeling and sensitivity analysis of lumped-parameter model of

[6] Saltelli A, Tarantola S, Chan KPS. A quantitative model-independent method for global

[7] Morris MD. Factorial sampling plans for preliminary computational experiments. Techno-

[8] Sobol I. Sensitivity estimates for nonlinear mathematical models. Matematicheskoe

[9] Homma T, Saltelli A. Importance measures in global sensitivity analysis of nonlinear

[10] Chen P, Quarteroni A, Rozza G. Simulation-based uncertainty quantification of human arterial network hemodynamics. International Journal for Numerical Methods in Bio-

[11] Gul R, Bernhard S. Optimal measurement locations for diagnosis of aortic stenoses and aneurysms in a lumped-parameter model of the systemic circulation using sensitivity

the human cardiovascular system. PhD Thesis. Germany: FU Berlin; 2016

sensitivity analysis of model output. Technometrics. 1999;41(1):39-56

Modelirovanie. 1990;2:112-118. In Russian, translated in English

models. Reliability Engineering and System Safety. 1996;52:117

analysis. International Journal of Biomathematics. 2017;10(8):1750116

Raheem Gul1

References

**Chapter 5**

**Provisional chapter**

**Biological Hypercomputation and Degrees of Freedom**

This chapter presents the idea of biological hypercomputation (BH) and discusses how and why it entails degrees of freedom. Crossing a biological and computational point of view, the claim is made that living beings cannot be considered as machines in any sense of the word, the arguments are provided, and the consequence is drawn: the complexity of life is the very process by which living beings gain degrees of freedom. The leading thread for the analysis here is the relationship between matter, energy and information.

Biological hypercomputation (BH) is the title that summarizes the fact that living systems are not machines, and hence, do not process information like any machine, no matter what. The leading thread, thereafter, is to understanding the very way in which living systems process information. For a living being processing information is after all a matter of life or death.

A normal understanding of the issue would claim that living systems receive information from the environment and then process it, so that they would be able to learn and adapt to the changing environment. This text claims that it is a wrong take. Instead, it will be argued that there is no information before the very act of processing, and as a consequence, there is no information after the processing. Reality is the outcome of the very process of information

To be sure, ranging from the various Turing Machines (TM)—u-TM, o-TM, d-TM, and many others, up to the various levels of life, say, from bacteria to living cells to plants to animals, etc., there is a growing process by which new and previously unforeseen degrees of freedom are reached. Reality, nature and the world can be seen as an increasing process of degrees of

**Biological Hypercomputation and Degrees of Freedom**

DOI: 10.5772/intechopen.73179

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution,

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

distribution, and reproduction in any medium, provided the original work is properly cited.

and reproduction in any medium, provided the original work is properly cited.

Carlos Eduardo Maldonado

**Abstract**

**1. Introduction**

processing.

Carlos Eduardo Maldonado

http://dx.doi.org/10.5772/intechopen.73179

Additional information is available at the end of the chapter

**Keywords:** life, computation, freedom, nature, complexity

Additional information is available at the end of the chapter

**Provisional chapter**

## **Biological Hypercomputation and Degrees of Freedom**

**Biological Hypercomputation and Degrees of Freedom**

DOI: 10.5772/intechopen.73179

#### Carlos Eduardo Maldonado Carlos Eduardo Maldonado Additional information is available at the end of the chapter

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/intechopen.73179

#### **Abstract**

This chapter presents the idea of biological hypercomputation (BH) and discusses how and why it entails degrees of freedom. Crossing a biological and computational point of view, the claim is made that living beings cannot be considered as machines in any sense of the word, the arguments are provided, and the consequence is drawn: the complexity of life is the very process by which living beings gain degrees of freedom. The leading thread for the analysis here is the relationship between matter, energy and information.

**Keywords:** life, computation, freedom, nature, complexity

#### **1. Introduction**

Biological hypercomputation (BH) is the title that summarizes the fact that living systems are not machines, and hence, do not process information like any machine, no matter what. The leading thread, thereafter, is to understanding the very way in which living systems process information. For a living being processing information is after all a matter of life or death.

A normal understanding of the issue would claim that living systems receive information from the environment and then process it, so that they would be able to learn and adapt to the changing environment. This text claims that it is a wrong take. Instead, it will be argued that there is no information before the very act of processing, and as a consequence, there is no information after the processing. Reality is the outcome of the very process of information processing.


**Dimension 4: Transformation**

Transformation: set concrete learning objectives and future action plans and approaches to initiate a new reflective cycle. Argumentation of these changes or the need for them.

**Table 3.** Elements of the rubric for students' identity construction assessment.


**2.3 The pharmacological activities of bacterial-derived secondary metabolites**

Pharmaceutical agents of bacterial origin include antibiotics, immunomodulators, nematicides, antitumor agents, coccidiostatic agents, enzyme inhibitors, and insecticides. Interestingly, *Escherichia coli* is used as a host in molecular biology for synthesis of recombinant proteins [76]. Selman Abraham Waksman, the father of antibiotics, discovered actinomycin; this effort was followed by the discovery of streptomycin in 1944 [77]. Other clinically important antibiotics derived from bacteria include bacitracin [78], polymyxin B [79], gentamicin [80], amphotericin b [81], tetracycline [82], erythromycin [83], rifamycin [84], vancomycin [82], neomycin [85], streptomycin [86], and chloramphenicol [87]. Etc.