# Modelling the Formation of Liver Zones within the Scope of Fractional Order Derivative

Authors: Abdon Atangana; Suares Clovis Oukouomi Noutchie
Journal: BioMed Research International
 (2014)

---

## Abstract

We develop and extend earlier results related to mathematical modelling of the liver formation zone by the adoption of noninteger order derivative. The hidden uncertainties in the model are captured and controlled thanks to the Caputo derivative. The stationary states are investigated and the time-dependent solution is approximated using two recent iteration methods. In particular, we discuss the convergence of these methods by constructing a suitable Hilbert space.


---

## Body

## 1. Introduction

One of the most important parts of any animal or human being body is the liver. The liver is a very important organ of vertebrates and many other flora and fauna [1]. Its properties include counting detoxification, protein synthesis, and production of biochemicals indispensable for absorption [2]. The liver is indispensable for continued existence; there is at this time no way to pay damages for the nonexistence of liver meaning in the long term, even though novel liver dialysis practices can be used in the temporary [1, 2]. This gland shows business of a main role in metabolism and has a number of functions in the body, including glycogen storage, decomposition of red blood cells, plasma protein synthesis, hormone production, and detoxification [3]. It is found below the diaphragm in the abdominal-pelvic area of the abdomen. It brings into being bile, an alkaline compound which aids in digestion via the emulsification of lipids. The liver’s extremely dedicated tissues normalize an extensive diversity of high-volume biochemical reactions, including the synthesis and breakdown of small and complex molecules, countless of which are essential for ordinary fundamental functions [3].The human liver is typically separated into two lobes (left and right), if viewed from the parietal surface, but if observed on the visceral surface it is divided into four lobes with the addition of the caudate and quadrate lobe [3]. Other anatomical landmarks exist, such as the ligamentum venosum (ligamentum of Arancio) and the round ligament (ligamentum Teres) that further divide the left side of the liver in two sections [4]. The falciform ligament is visible on the front (anterior side) of the liver. This divides the liver into a left anatomical lobe and a right anatomical lobe [1–4]. Two most important types of cells inhabit the liver lobes: karat parenchymal and nonparenchymal cells. 80% of the liver volume is engaged by parenchymal cells normally referred to as hepatocytes. Nonparenchymal cells comprise 40% of the total number of liver cells but only 6.5% of its volume [3]. Sinusoidal hepatic endothelial cells, Kupffer cells, and hepatic stellate cells are some of the nonparenchymal cells that line the liver sinusoid [3]. A schematic representation of anatomy of the biliary tree, liver, and gall bladder is shown in Figure 1.Figure 1
Anatomy of the biliary tree, liver, and gall bladder (From Jiju Kurian Punnoose, 10/12/2007).More than a few metabolic occupations of the liver have been established to be prearranged in spatial zones arranged in relation to the direction of hepatic blood flow, in such a way that a quantity of enzymes operate approximately together upstream others. A group of the most eminent scholars who studied this phenomenon are Bass et al. [5]. In one of their investigations, they qualified such distributions of enzymes activities to distributions of cell types [5]. For the simplest case of two enzymes, there are two corresponding cell types, each containing only one of the enzymes; separate metabolic zones occur when all cells of one type are located upstream all cells of the other type. Additionally, it was stated in [5] that each cell type reproduces itself by division.The mathematical model of the formation of liver zone was discussed in [5]. However in order to accommodate readers that are not in the field of mathematical medicine, we will discuss for expediency the main steps in its derivation in the next section.In the recent decades it was revealed by several scholars dealing with real world problems that the models using the concept of noninteger order derivatives are more suitable in prediction than those with ordinary derivatives. One of the goals of this paper is the analysis of the model within the folder of fractional calculus including the steady state analysis, a possible analytical solution using the recent development of analytical methods.


## 2. Mathematical Formulation

The mathematical model describing the formation of liver zones is a system of nonlinear integropartial differential equations. Approximately 1100 mL of blood flows from the portal vein into the liver sinusoids each minute, and approximately an additional 350 mL flows into the sinusoids from the hepatic artery, and the total averaging is about 1450 mL/min. This amounts to about 29% of the resting cardiac output [6]. As the many capillaries comprising the liver are similar and act essentially in parallel, Bass et al. [5] modelled a representative capillary lined with cells of two kinds. In the investigation done by Bass et al, it was recommended to put the x-axis along the blood flow, with bay at x
=
0 and outlet at x
=
L. The density of cells of the first and second kinds is defined by ρ
1
(
x
,
t
) and ρ
2
(
x
,
t
), respectively, as a continuous representation of the number of cells of the first kind and second per unit length of capillary at time t at the position x. It was suggested in their investigation that the total cell density ρ
1
+
ρ
2 cannot go beyond a fixed highest density σ of cell sites, as division of the cell is limited by the well-known phenomenon of contact inhibition. The mathematical formulation describing the above phenomena is given as follows:
(1)
∂
ρ
1
∂
t
=
K
1
ρ
1
(
σ
-
ρ
1
-
ρ
2
)
-
β
1
(
c
)
ρ
1
,
∂
ρ
2
∂
t
=
K
2
ρ
2
(
σ
-
ρ
1
-
ρ
2
)
-
β
2
(
c
)
ρ
2
,
c
=
c
0
-
1
f
∫
0
x
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
.
Here k
1 and k
2 are positive constants and c is the concentration of a controlling blood-borne substance. The above system was reduced in [5] as follows:
(2)
∂
ρ
1
∂
t
=
K
1
ρ
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
ρ
1
γ
1
f
∫
0
x
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
,
∂
ρ
2
∂
t
=
K
2
ρ
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
ρ
2
γ
2
f
∫
0
x
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
.
In this paper we are not going to consider the above system since in recent decades it was revealed by several scholars dealing with real world problems that the models using the concept of noninteger order derivatives are more suitable in prediction than those with ordinary derivative [7–10]. In order to take into account all uncertainties that can be associated with the formation of the liver zone in human or animal body, we will convert the ordinary derivative to the concept of noninteger order derivative as follows:
(3)
∂
α
ρ
1
∂
t
α
=
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
1
,
∂
α
ρ
2
∂
t
α
=
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
2
,
where 0
<
α
≤
1, μ
i
i
=
1,2, are death rate of kinds 1 and 2, respectively, and γ
i
i
=
1,2 are constants. This conversion also implies that the global concentration of a controlling blood-borne substance is presented as
(4)
c
(
x
,
t
,
α
)
=
c
0
-
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
,
where Γ is the gamma function defined as
(5)
Γ
(
x
)
=
∫
0
∞
t
x
-
1
e
-
t
d
t
.
It is worth noting that if α
=
1, we revert the conventional model of liver zone formation. The fractional derivative used here is in Caputo sense. There exists more than one definition of fractional derivative in the literature, and we will list the most common ones. The most popular ones are the Riemann-Liouville and the Caputo derivatives [8–10]. For Caputo derivative we have
(6)
D
0
C
x
α
(
f
(
x
)
)
=
1
Γ
(
n
-
α
)
∫
0
x
(
x
-
t
)
n
-
α
-
1
d
n
f
(
t
)
d
t
n
d
t
.
For the case of Riemann-Liouville we have the following definition:
(7)
D
x
α
(
f
(
x
)
)
=
1
Γ
(
n
-
α
)
d
n
d
x
n
∫
0
x
(
x
-
t
)
n
-
α
-
1
f
(
t
)
d
t
.
It was shown in [7–9] that the Caputo type is suitable for modelling real world problem, and based on this we will use throughout this paper the Caputo derivative.


## 3. Steady State Analysis

This section is devoted to the discussions concerning the stationary solutions, which means that we will consider that the densities are not depending on time. One of the enjoyable properties of the Caputo fractional derivative is that the fractional derivative of a constant is zero. Now based on that property, we assume that the densities of first and second kind do not depend on time; (3) will be reduced to
(8)
0
=
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
+
k
2
ρ
2
*
]
K
1
(
σ
-
ρ
1
*
-
ρ
2
*
)
-
μ
1
-
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
+
k
2
ρ
2
*
]
d
φ
]
ρ
1
,
0
=
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
+
k
2
ρ
2
*
]
K
2
(
σ
-
ρ
1
*
-
ρ
2
*
)
-
μ
2
-
γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
+
k
2
ρ
2
*
]
d
φ
]
ρ
2
.
This simply implies
(9)
K
1
(
σ
-
ρ
1
*
-
ρ
2
*
)
-
μ
1
-

γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*

+

k
2
ρ
2
*
]
d
φ
with
ρ
2
=
0
,
(10a)
K
2
(
σ
-
ρ
1
*
-
ρ
2
*
)
-
μ
2
-

γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*

+

k
2
ρ
2
*
]
d
φ
,
with
ρ
1
=
0
,
(10b)
ρ
1
-
ρ
2
=
0
,
(11)
K
1
(
σ
-
ρ
1
*
-
ρ
2
*
)
-
μ
1
-

γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*

+

k
2
ρ
2
*
]
d
φ
with
ρ
2
*


≠

0
,
K
2
(
σ
-
ρ
1
*
-
ρ
2
*
)
-
μ
2
-
γ
2
f
Γ
(
α
)
∫
0
x

(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*

+

k
2
ρ
2
*
]
d
φ
,
with
ρ
1
*

≠

0
.
We will start with the stationary solution of system (9). We will employ two different analytical techniques to get to the bottom of the system. We will start with the Laplace transform method:
(12)
K
1
(
σ
-
ρ
1
*
)
-
μ
1
-
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
]
d
φ
=
0
.
Now if we put for simplicity d
1
=
σ
-
(
μ
1
/
K
1
)
,
a
1
=
k
1
γ
1
/
f
K
1 with rearranging we obtain the following:
(13)
ρ
1
*
=
d
1
-
a
1
1
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
ρ
1
*
]
d
φ
.
We will first show the properties of the Laplace transform in fractional calculus [8–10].The Laplace transform of the functionf is defined as follows:
(14)
L
(
f
(
x
)
)
(
s
)
=
∫
0
∞
e
-
s
x
f
(
x
)
d
x
.
Let us observe the Laplace transform of the fractional derivative with Caputo
(15)
L
(
D
0
C
t
α
f
(
x
)
)
(
s
)
=
s
α
F
(
s
)
-
∑
k
=
0
n
-
1
s
α
-
k
-
1
f
(
k
)
(
0
)
,
(
n
-
1
<
α
≤
n
)
.
The above uses the usual initial conditions or values of the functions. Also the Laplace transform of fractional integral is given as
(16)
L
(
J
t
α
f
(
x
)
)
(
s
)
=
s
-
α
F
(
s
)
,
with
J
α
f
(
x
)
=
1
Γ
(
α
)
∫
0
x
(
x
-
t
)
α
-
1
f
(
t
)
d
t
,
α
>
0
,
x
>
0
.
F is the Laplace transform of f and s the Laplace variable. Therefore applying the Laplace transform on both sides of (12), we arrive at
(17)
L
{
ρ
1
*
(
x
)
}
(
s
)
=
d
1
s
+
a
1
s
1
-
α
.
Now applying the inverse-Laplace transform operator on both sides of (17) we arrive at
(18)
ρ
1
*
(
x
)
=
L
-
1
[
d
1
s
+
a
1
s
1
-
α
]
(
x
)
=
d
1
E
α
[
-
(
a
1
1
/
α
x
)
α
]
,
where E
α is the so-called Mittag-Leffler function defined as
(19)
E
α
(
x
)
=
∑
n
=
0
∞
x
n
Γ
(
α
n
+
1
)
.
The above solution can be retrieved by using new recent development of iteration methods. This method is being proposed to avoid the stability of Laplace transform method on one hand. Also the Laplace transform method is only useful in the case of linear equation; however, the second method is used for both linear and nonlinear equations.In this method we assume that the solution of (13) is in the form of
(20)
ρ
1
*
(
x
)
=
lim
⁡
p
→
1
⁡
∑
n
=
0
∞
p
n
ρ
1
n
(
x
)
.
Now substituting the above in (13) we obtain
(21)
∑
n
=
0
∞
p
n
ρ
1
n
(
x
)
=
d
1
-
a
1
p
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
∑
n
=
0
∞
p
n
ρ
1
n
(
x
)
]
d
φ
.
Comparing terms of the same power of p yields the following iteration:
(22)
ρ
10
=
d
1
,
ρ
11
(
x
)
=
-
a
1
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
ρ
10
(
φ
)
]
d
φ
,
ρ
12
(
x
)
=
-
a
1
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
ρ
11
(
φ
)
]
d
φ
,
ρ
1
n
(
x
)
=
-
a
1
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
ρ
1
(
n
-
1
)
(
φ
)
]
d
φ
.
Integrating the above set of integral equations, we arrive at the following result:
(23)
ρ
10
=
d
1
,
ρ
11
(
x
)
=
-
d
1
a
1
x
α
Γ
(
α
+
1
)
,
ρ
12
(
x
)
=
d
1
(
-
a
1
)
2
(
x
)
2
α
Γ
(
2
α
+
1
)
,
ρ
1
n
(
x
)
=
d
1
(
-
a
1
)
n
(
x
)
n
α
Γ
(
n
α
+
1
)
.
Implying the stationary solution of density of the first kind is given as
(24)
ρ
1
*
(
x
)
=
d
1
∑
n
=
0
∞
(
-
a
1
)
n
(
x
)
n
α
Γ
(
n
α
+
1
)
=
d
1
E
α
[
-
(
(
a
1
)
1
/
α
x
)
α
]
.
It was argued in the work done in [5] that there exists a certain point x
=
x
* above which the density of the second kind ρ
2
(
x
) will approach a stationary, which can be determined by solving (25). Therefore in that interval we have the following equation:
(25)
K
2
(
σ
-
ρ
2
*
)
-
μ
2
-
γ
2
f
Γ
(
α
)
∫
0
x
*
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
]
d
φ
-
γ
2
f
Γ
(
α
)
∫
x
*
x
(
x
-
φ
)
α
-
1
[
k
2
ρ
2
*
]
d
φ
=
0
.
We can rearrange the above equation as follows:
(26)
ρ
2
*
(
x
)
=
F
(
x
*
)
-
γ
2
k
2
f
K
2
Γ
(
α
)
∫
x
*
x
(
x
-
φ
)
α
-
1
[
ρ
2
*
]
d
φ
,
F
(
x
*
)
=
σ
-
μ
2
K
2
-
γ
2
f
Γ
(
α
)
K
2
∫
0
x
*
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
]
d
φ
.
For simplicity, let us put
(27)
γ
2
k
2
f
K
2
=
a
2
,
ρ
2
*
(
x
)
=
F
(
x
*
)
-
a
2
Γ
(
α
)
∫
x
*
x
(
x
-
φ
)
α
-
1
[
ρ
2
*
]
d
φ
.
Using a similar method, we obtain
(28)
ρ
20
=
F
(
x
*
)
,
ρ
21
(
x
)
=
-
a
2
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
ρ
20
(
φ
)
]
d
φ
,
ρ
22
(
x
)
=
-
a
2
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
ρ
21
(
φ
)
]
d
φ
,
ρ
2
n
(
x
)
=
-
a
2
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
ρ
2
(
n
-
1
)
(
φ
)
]
d
φ
.
Integrating the above set of integral equations, we arrive at the following result:
(29)
ρ
20
=
F
(
x
*
)
ρ
21
(
x
)
=
-
F
(
x
*
)
a
2
(
x
-
x
*
)
α
Γ
(
α
+
1
)
,
ρ
22
(
x
)
=
F
(
x
*
)
(
-
a
2
)
2
(
x
-
x
*
)
2
α
Γ
(
2
α
+
1
)
,
ρ
2
n
(
x
)
=
F
(
x
*
)
(
-
a
2
)
n
(
x
-
x
*
)
n
α
Γ
(
n
α
+
1
)
.
Implying the stationary solution of density of the second kind is given as
(30)
ρ
1
*
(
x
)
=
F
(
x
*
)
∑
n
=
0
∞
(
-
a
2
)
n
(
x
-
x
*
)
n
α
Γ
(
n
α
+
1
)
=
F
(
x
*
)
E
α
[
-
(
(
a
2
)
1
/
α
(
x
-
x
*
)
)
α
]
.
Our next concern is to provide the stationary solutions for the following system:
(31)
K
1
(
σ
-
ρ
1
*
-
ρ
2
*
)
-
μ
1
-
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
+
k
2
ρ
2
*
]
d
φ
=
0
with
ρ
2
*
≠
0
,
K
2
(
σ
-
ρ
1
*
-
ρ
2
*
)
-
μ
2
-
γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
+
k
2
ρ
2
*
]
d
φ
=
0
,
with
ρ
1
*
≠
0
.
Dividing the first equation of the system by K
1 and the second by K
2 and subtracting the result of the first from the result of the second, we obtain
(32)
1
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
*
+
k
2
ρ
2
*
]
d
φ
=
(
μ
1
/
K
1
)
-
(
μ
2
/
K
2
)
(
γ
2
/
f
K
2
)
-
(
γ
1
/
f
K
1
)
.
Now replacing the above result in the system, we arrive at
(33)
σ
-
ρ
1
*
-
ρ
2
*
=
μ
1
+
γ
1
f
K
1
(
μ
1
/
K
1
)
-
(
μ
2
/
K
2
)
(
γ
2
/
f
K
2
)
-
(
γ
1
/
f
K
1
)
=
μ
2
+
γ
2
f
K
2
(
μ
1
/
K
1
)
-
(
μ
2
/
K
2
)
(
γ
2
/
f
K
2
)
-
(
γ
1
/
f
K
1
)
.
A condition for the above to hold is that (
γ
2
/
f
K
2
)
-
(
γ
1
/
f
K
1
)
≠
0
,
(
μ
1
/
K
1
)
-
(
μ
2
/
K
2
) must have the same sign with (
γ
2
/
f
K
2
)
-
(
γ
1
/
f
K
1
), and the ultimate best is that
(34)
μ
1
K
1
>
μ
2
K
2
,
γ
2
f
K
2
>
γ
1
f
K
1
.
Thus, it is confirmed in this analysis that, even with stationary solutions,
(35)
σ
-
ρ
1
*
-
ρ
2
*
>
0
.
We will now continue our investigation by presenting methods to find the solutions of system (3). This is therefore done in the next section.


## 4. Finding Approximate Global Solution

Many real world problems have been modelled via differential equations. It is important to point out that once these physical problems are converted into mathematical equation, their solutions are further used to predict the future behaviors of these physical problems. It is therefore very important to develop analytical methods in order to find the solutions. Some of these equations are very difficult to handle analytically due to their complexity especially in the case of fractional calculus. In the last years, many scientists focused their regards on setting up techniques that can be used to solve these equations. But for our purpose, we will use only two of these techniques to present an approximate solution of (3).


### 4.1. Iteration Method Using Laplace Transform

This method was proposed by Khan and Wu to avoid calculation of Adomian polynomial and the use of correction function of the variational iteration method; the reference of this work is in [11]. One can also find the methodology of this efficient method in there, but here we will use it only. Now applying the Laplace transform on both sides of the below system, we arrive at
(36)
∂
α
ρ
1
∂
t
α
=
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
1
,
∂
α
ρ
2
∂
t
α
=
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
2
and the following
(37)
s
α
L
[
ρ
1
(
x
,
t
)
]
(
s
)
=
L
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
1
,
s
α
L
[
ρ
2
(
x
,
t
)
]
(
s
)
=
L
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1

+

k
2
ρ
2
]
d
φ
]
ρ
2
.
For simplicity, let us put
(38)
L
[
ρ
1,2
(
x
,
t
)
]
(
s
)
=
ρ
1,2
ˇ
.
Now taking the inverse Laplace operator on both sides of (37), we arrive at
(39)
ρ
1
(
x
,
t
)
=
ρ
1
(
x
,
0
)
+
L
-
1
[
1
s
α
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
γ
1
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
1
]
,
ρ
2
(
x
,
t
)
=
ρ
2
(
x
,
0
)
+
L
-
1
[
1
s
α
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
γ
2
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
2
]
.
Since we are not yet sure of the exact solution of the above system, we will assume that the solution can be represented in series form according to the Poincare idea; that is,
(40)
ρ
1
(
x
,
t
)
=
∑
n
=
0
∞
p
n
ρ
1
n
(
x
,
t
)
,
ρ
2
(
x
,
t
)
=
∑
n
=
0
∞
p
n
ρ
2
n
(
x
,
t
)
.
The above can be substituted in (39); however, comparing terms of the same power of p we arrived at the following set:
(41)
ρ
10
(
x
,
t
)
=
ρ
1
(
x
,
0
)
,
ρ
20
(
x
,
t
)
=
ρ
2
(
x
,
0
)
,
ρ
11
(
x
,
t
)
=
+
L
-
1
[
1
s
α
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
K
1
(
σ
-
ρ
10
-
ρ
20
)
-
μ
1
-
γ
1
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
d
φ
]
ρ
10
]
,
ρ
21
(
x
,
t
)
=
+
L
-
1
[
1
s
α
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
K
1
(
σ
-
ρ
10
-
ρ
20
)
-
μ
2
-
γ
2
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
d
φ
]
ρ
20
]
.
In general for n greater than 1 we can obtain the rest of the terms with the following recursive formula:
(42)
ρ
1
,
n
+
1
(
x
,
t
)
=
L
-
1
[
1
s
α
L
[
K
1
σ
ρ
1
,
n
-
K
1
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
1
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j

-

μ
1
ρ
1
,
n

-

ρ
1
n
J
1
α
n
(
x
,
t
)
]
]
,
ρ
2
,
n
+
1
(
x
,
t
)
=
L
-
1
[
1
s
α
L
[
K
2
σ
ρ
2
,
n
-
K
2
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
2
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
2
α
n
(
x
,
t
)
]
]
with of course
(43)
J
2
α
n
(
x
,
t
)
=
γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
n
+
k
2
ρ
2
n
]
d
φ
,
J
1
α
n
(
x
,
t
)
=
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
n
+
k
2
ρ
2
n
]
d
φ
.
We will summarize the entire procedure in the following algorithm.Algorithm 1.
(i) Input:I
1
(
x
)
,
I
2
(
x
)—as initial conditions, meaning in our case ρ
1
(
x
,
0
) and ρ
1
(
x
,
0
), j—number terms in the rough calculation.
(ii) Output:ρ
1
approx
(
x
,
t
)
,
ρ
2
approx
(
x
,
t
) the approximate solutions.Step 1.
Putρ
10
(
x
,
t
)
=
I
1
(
x
)
,
ρ
20
(
x
,
t
)
=
I
2
(
x
) and ρ
1
approx
(
x
,
t
)
=
ρ
10
(
x
,
t
)
,
ρ
2
approx
(
x
,
t
)
=
I
2
(
x
).Step 2.
Forj
=
0 to n
-
1 do Step 3, Step 4, and Step 5.Step 3.
Compute(44)
M
1
n
=
K
1
σ
ρ
1
,
n
-
K
1
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
1
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
1
α
n
(
x
,
t
)
,
M
2
n
=
K
2
σ
ρ
2
,
n
-
K
2
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
2
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
2
α
n
(
x
,
t
)
.Step 4.
Compute(45)
ρ
1
,
n
+
1
(
x
,
t
)
=
L
-
1
[
1
s
α
L
[
M
1
n
]
]
,
ρ
2
,
n
+
1
(
x
,
t
)
=
L
-
1
[
1
s
α
L
[
M
2
n
]
]
.Step 5.
Computeρ
1
approx
(
x
,
t
)
=
ρ
1
approx
(
x
,
t
)
+
ρ
1
,
n
+
1
(
x
,
t
), ρ
2
approx
(
x
,
t
)
=
ρ
2
approx
(
x
,
t
)
+
ρ
2
,
n
+
1
(
x
,
t
).
Stop.


### 4.2. Iteration Method Using the Integral Transform

This version of homotopy was first proposed by Abdon Atangana and used to solve the groundwater flow equation in the confined aquifer [12, 13]. This version uses the idea of integral transform coupling with the Poincare series that was previously proposed by He [14]. This version is easier to use as one can see in this case.The first step in this approach is to apply the inverse operator of∂
α
/
∂
t
α on both sides of equation to obtain the following equation:
(46)
ρ
1
=
ρ
1
(
x
,
0
)
+
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
γ
1
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
1
d
φ
,
ρ
2
=
ρ
2
(
x
,
0
)
+
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
γ
2
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
2
d
φ
.
Since it is not certain that the above has exact solution, the second step of the method proposes that the solution should be in the series of Poincare meaning:
(47)
ρ
1
(
x
,
t
)
=
∑
n
=
0
∞
p
n
ρ
1
n
(
x
,
t
)
,
ρ
2
(
x
,
t
)
=
∑
n
=
0
∞
p
n
ρ
2
n
(
x
,
t
)
.
Then ρ
1,2 can be included in (46) and after including the imbedding parameter p and comparing terms of the same power of p we arrive at the following set of equations:
(48)
ρ
10
(
x
,
t
)
=
ρ
1
(
x
,
0
)
,
ρ
20
(
x
,
t
)
=
ρ
2
(
x
,
0
)
,
ρ
11
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
K
1
(
σ
-
ρ
10
-
ρ
20
)
-
μ
1
-
γ
1
f
Γ
(
α
)
×

∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10

+

k
2
ρ
20
]
d
φ
]
ρ
10
d
ψ
,
ρ
21
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
K
1
(
σ
-
ρ
10
-
ρ
20
)
-
μ
2
-
γ
2
f
Γ
(
α
)
×

∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10

+

k
2
ρ
20
]
d
φ
]
ρ
20
d
ψ
.
In general for n greater than 1 we can obtain the rest of the terms with the following iteration formula:
(49)
ρ
1
,
n
+
1
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
K
1
σ
ρ
1
,
n
-
K
1
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
1
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
1
α
n
(
x
,
t
)
]
d
ψ
,
ρ
2
,
n
+
1
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
K
2
σ
ρ
2
,
n
-
K
2
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
2
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
2
α
n
(
x
,
φ
)
]
d
ψ
.
As soon as one is provided with initial condition, this algorithm can be implemented easily. To make this simple, we provide the summary of the procedure for this case in the following algorithm.Algorithm 2.
(i) Put in:I
1
(
x
)
,
I
2
(
x
)—as preliminary conditions, meaning in our case ρ
1
(
x
,
0
) and ρ
1
(
x
,
0
), k—number terms in the rough calculation.
(ii) Output:ρ
1
approx
(
x
,
t
)
,
ρ
2
approx
(
x
,
t
) the approximate solutions.Step 1.
Putρ
10
(
x
,
t
)
=
I
1
(
x
)
,
ρ
20
(
x
,
t
)
=
I
2
(
x
) and ρ
1
approx
(
x
,
t
)
=
ρ
10
(
x
,
t
)
,
ρ
2
approx
(
x
,
t
)
=
I
2
(
x
).Step 2.
Fork
=
0 to n
-
1 do Step 3, Step 4, and Step 5.Step 3.
Compute(50)
L
1
n
=
K
1
σ
ρ
1
,
n
-
K
1
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
1
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
1
α
n
(
x
,
t
)
,
L
2
n
=
K
2
σ
ρ
2
,
n
-
K
2
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
2
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
2
α
n
(
x
,
t
)
.Step 4.
Compute(51)
ρ
1
,
n
+
1
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
L
1
n
d
ψ
,
ρ
2
,
n
+
1
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
L
2
n
d
ψ
.Step 5.
Computeρ
1
approx
(
x
,
t
)
=
ρ
1
approx
(
x
,
t
)
+
ρ
1
,
n
+
1
(
x
,
t
), ρ
2
approx
(
x
,
t
)
=
ρ
2
approx
(
x
,
t
)
+
ρ
2
,
n
+
1
(
x
,
t
).
Stop.Remark 3.
A direct comparison of both methods shows that they all use integral transform and the Poincare series. However in the first method, one will first apply the Laplace transform and again apply the inverse-Laplace operator which sometime can be very difficult to handle due to the instability of the Laplace transform operator. But the second method uses just a simple integral transform that can be sometime implemented with the new powerful software. We therefore conclude that the second method is much easier to implement rather than the first method. Nevertheless both methods have been proven to be efficient.In the recent decade it was observed by several editors that many scholars using iteration method to find approximate or exact solutions of complicated equations were not investigating the stability and the convergence of their iterations. This observation led us to think that this field was much easier, and other editors said it was a child play. In order to avoid classifying this paper under this class, we will present without loss of generality the convergence of our iteration number 2 and this is done in the next section.


### 4.3. Convergence Investigation

In this section, we examine the convergence property of the approximated solution for the formation of liver zone equation. Let us consider the formation of liver zone equation in the Hilbert spaceH
=
L
2
(
(
η
,
λ
)
×
[
0
,
T
]
) defined as
(52)
H

=


{
(
u
,
v
)
:
(
η
,
λ
)

×

[
0
,
T
]
,
with
∫
(
t
-
φ
)
α
-
1
Γ
(
α
)
u
v
d
ι
d
κ

<

∞
}
.
Then the operator is of the form
(53)
H
(
u
,
v
)
=
u
(
a
(
σ
-
u
-
v
)
-
μ
-
b
∫
0
x
[
c
u
+
d
v
]
d
φ
)
.
That can be reduced to the following operator due to the physical properties of the problem under investigation:
(54)
H
(
u
)
=
u
(
a
(
σ
-
u
-
u
)
-
μ
-
b
∫
0
x
[
c
u
+
d
u
]
d
φ
)
.
The homotopy decomposition method is convergent if the following conditions are satisfiedHypothesis 1.
It is possible for us to find a positive constant, say,F, such that the inner product holds in the Hilbert space H
(55)
(
H
(
u
)
-
H
(
v
)
,
u
-
v
)
≥
F
∥
u
-
v
∥
,
∀
v
,
u
∈
H
.Hypothesis 2.
As far as for allv
,
u
∈
H are bounded this implies, we can find a positive constant, say, C, such that ∥
u
∥
,
∥
v
∥
≤
C, then we can find Φ
(
C
)
>
0 such that
(56)
(
H
(
u
)
-
H
(
v
)
,
g
)
≥
Φ
(
C
)
∥
u
-
v
∥
∥
z
∥
,
∀
z
∈
H
.
We can therefore state the following theorem for the sufficient condition for the convergence of the generalized equation.Theorem 4.
Let us consider(57)
H
(
u
)
=
∂
α
u
∂
t
α
=
u
(
a
(
σ
-
u
-
u
)
-
μ
-
b
∫
0
x
[
c
u
+
d
u
]
d
φ
)
,
w
i
t
h
a
σ
>
2
a
+
b
2
(
c
+
d
)
and consider the initial and boundary condition for the generalized equation; then, the homotopy decomposition method leads to a special solution of system (4).We will present the proof of this theorem by just verifying Hypothesises1 and 2.Proof.
We have that(58)
H
(
u
)
-
H
(
v
)
=
a
σ
(
u
-
v
)
-
2
a
(
u
2
-
v
2
)
-
μ
(
v
-
u
)
-
b
1
u
∫
0
x
[
u
]
d
φ
+
b
1
v
∫
0
x
[
v
]
d
φ
.
The above equation can be roughly converted to
(59)
H
(
u
)
-
H
(
v
)
=
a
σ
(
u
-
v
)
-
2
a
(
u
2
-
v
2
)
-
μ
(
u
-
v
)
-
b
1
2
d
(
u
2
-
v
2
)
d
x
=
a
σ
(
u
-
v
)
-
2
a
(
u
-
v
)
(
u
+
v
)
-
b
1
2
d
[
(
u
-
v
)
(
u
+
v
)
]
d
x
.
However according to the physical property of this problem, we have that ρ
1
+
ρ
2 the total cell density cannot exceed some fixed maximum density σ of cell sites, as limited by familiar phenomenon of contact inhibition; thus, we are urged to conclude that we can find a positive constant, say, h, such that u
+
v
=
h, and then,
(60)
=
a
σ
(
u
-
v
)
-
2
a
(
u
-
v
)
(
u
+
v
)
-
b
1
2
d
[
(
u
-
v
)
(
u
+
v
)
]
d
x
can be converted to
(61)
=
a
σ
(
u
-
v
)
-
2
a
(
u
-
v
)
h
-
b
1
2
d
[
(
u
-
v
)
h
]
d
x
.
With the above in hand, we have the following inner product:
(62)
(
H
(
u
)
-
H
(
v
)
,
u
-
v
)
=
a
σ
(
(
u
-
v
)
,
u
-
v
)
-
2
a
(
(
u
-
v
)
,
u
-
v
)
h
-
h
b
1
2
d
[
(
(
u
-
v
)
,
u
-
v
)
]
d
x
.
According to the physical property of this problem, u
,
v are bounded; therefore, we can find a positive constant F such that (
u
,
u
)
,
(
v
,
v
)
<
F
2. It follows by the use of Schwartz inequality that
(63)
d
[
(
(
u
-
v
)
,
u
-
v
)
]
d
x
≤
∥
(
u
-
v
)
x
∥
∥
u
-
v
∥
.
Now since there exists a positive constant ω such that ∥
(
u
-
v
)
x
∥
≤
ω
∥
u
-
v
∥ then
(64)
d
[
(
(
u
-
v
)
,
u
-
v
)
]
d
x
≤
ω
∥
u
-
v
∥
2
.
On the other hand, we have that
(65)
(
(
u
-
v
)
,
u
-
v
)
=
∥
u
-
v
∥
2
.
Now making a direct substitution of (64) and (65) into (62) we arrive at
(66)
(
H
(
u
)
-
H
(
v
)
,
u
-
v
)
≥
F
∥
u
-
v
∥
2
with F
=
a
σ
-
2
a
h
-
(
h
b
1
/
2
). Therefore Hypothesis 1 is true for the generalized version. We will now present the proof for Hypothesis 2, as we are directly computing
(67)
(
H
(
u
)
-
H
(
v
)
,
z
)
=
a
σ
(
(
u
-
v
)
,
z
)
-

2
a
(
(
u
-
v
)
,
z
)
h

-
h
b
1
2
d
[
(
(
u
-
v
)
,
z
)
]
d
x
≥
Φ
(
C
)
∥
u
-
v
∥
∥
z
∥
.
Therefore we can see that Hypothesis 2 is correct as well. We can conclude without fear that the homotopy decomposition method works perfectly for this generalized equation.


## 4.1. Iteration Method Using Laplace Transform

This method was proposed by Khan and Wu to avoid calculation of Adomian polynomial and the use of correction function of the variational iteration method; the reference of this work is in [11]. One can also find the methodology of this efficient method in there, but here we will use it only. Now applying the Laplace transform on both sides of the below system, we arrive at
(36)
∂
α
ρ
1
∂
t
α
=
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
1
,
∂
α
ρ
2
∂
t
α
=
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
2
and the following
(37)
s
α
L
[
ρ
1
(
x
,
t
)
]
(
s
)
=
L
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
1
,
s
α
L
[
ρ
2
(
x
,
t
)
]
(
s
)
=
L
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1

+

k
2
ρ
2
]
d
φ
]
ρ
2
.
For simplicity, let us put
(38)
L
[
ρ
1,2
(
x
,
t
)
]
(
s
)
=
ρ
1,2
ˇ
.
Now taking the inverse Laplace operator on both sides of (37), we arrive at
(39)
ρ
1
(
x
,
t
)
=
ρ
1
(
x
,
0
)
+
L
-
1
[
1
s
α
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
γ
1
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
1
]
,
ρ
2
(
x
,
t
)
=
ρ
2
(
x
,
0
)
+
L
-
1
[
1
s
α
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
γ
2
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
2
]
.
Since we are not yet sure of the exact solution of the above system, we will assume that the solution can be represented in series form according to the Poincare idea; that is,
(40)
ρ
1
(
x
,
t
)
=
∑
n
=
0
∞
p
n
ρ
1
n
(
x
,
t
)
,
ρ
2
(
x
,
t
)
=
∑
n
=
0
∞
p
n
ρ
2
n
(
x
,
t
)
.
The above can be substituted in (39); however, comparing terms of the same power of p we arrived at the following set:
(41)
ρ
10
(
x
,
t
)
=
ρ
1
(
x
,
0
)
,
ρ
20
(
x
,
t
)
=
ρ
2
(
x
,
0
)
,
ρ
11
(
x
,
t
)
=
+
L
-
1
[
1
s
α
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
K
1
(
σ
-
ρ
10
-
ρ
20
)
-
μ
1
-
γ
1
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
d
φ
]
ρ
10
]
,
ρ
21
(
x
,
t
)
=
+
L
-
1
[
1
s
α
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
K
1
(
σ
-
ρ
10
-
ρ
20
)
-
μ
2
-
γ
2
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
d
φ
]
ρ
20
]
.
In general for n greater than 1 we can obtain the rest of the terms with the following recursive formula:
(42)
ρ
1
,
n
+
1
(
x
,
t
)
=
L
-
1
[
1
s
α
L
[
K
1
σ
ρ
1
,
n
-
K
1
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
1
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j

-

μ
1
ρ
1
,
n

-

ρ
1
n
J
1
α
n
(
x
,
t
)
]
]
,
ρ
2
,
n
+
1
(
x
,
t
)
=
L
-
1
[
1
s
α
L
[
K
2
σ
ρ
2
,
n
-
K
2
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
2
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
2
α
n
(
x
,
t
)
]
]
with of course
(43)
J
2
α
n
(
x
,
t
)
=
γ
2
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
n
+
k
2
ρ
2
n
]
d
φ
,
J
1
α
n
(
x
,
t
)
=
γ
1
f
Γ
(
α
)
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
n
+
k
2
ρ
2
n
]
d
φ
.
We will summarize the entire procedure in the following algorithm.Algorithm 1.
(i) Input:I
1
(
x
)
,
I
2
(
x
)—as initial conditions, meaning in our case ρ
1
(
x
,
0
) and ρ
1
(
x
,
0
), j—number terms in the rough calculation.
(ii) Output:ρ
1
approx
(
x
,
t
)
,
ρ
2
approx
(
x
,
t
) the approximate solutions.Step 1.
Putρ
10
(
x
,
t
)
=
I
1
(
x
)
,
ρ
20
(
x
,
t
)
=
I
2
(
x
) and ρ
1
approx
(
x
,
t
)
=
ρ
10
(
x
,
t
)
,
ρ
2
approx
(
x
,
t
)
=
I
2
(
x
).Step 2.
Forj
=
0 to n
-
1 do Step 3, Step 4, and Step 5.Step 3.
Compute(44)
M
1
n
=
K
1
σ
ρ
1
,
n
-
K
1
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
1
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
1
α
n
(
x
,
t
)
,
M
2
n
=
K
2
σ
ρ
2
,
n
-
K
2
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
2
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
2
α
n
(
x
,
t
)
.Step 4.
Compute(45)
ρ
1
,
n
+
1
(
x
,
t
)
=
L
-
1
[
1
s
α
L
[
M
1
n
]
]
,
ρ
2
,
n
+
1
(
x
,
t
)
=
L
-
1
[
1
s
α
L
[
M
2
n
]
]
.Step 5.
Computeρ
1
approx
(
x
,
t
)
=
ρ
1
approx
(
x
,
t
)
+
ρ
1
,
n
+
1
(
x
,
t
), ρ
2
approx
(
x
,
t
)
=
ρ
2
approx
(
x
,
t
)
+
ρ
2
,
n
+
1
(
x
,
t
).
Stop.


## 4.2. Iteration Method Using the Integral Transform

This version of homotopy was first proposed by Abdon Atangana and used to solve the groundwater flow equation in the confined aquifer [12, 13]. This version uses the idea of integral transform coupling with the Poincare series that was previously proposed by He [14]. This version is easier to use as one can see in this case.The first step in this approach is to apply the inverse operator of∂
α
/
∂
t
α on both sides of equation to obtain the following equation:
(46)
ρ
1
=
ρ
1
(
x
,
0
)
+
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
1
(
σ
-
ρ
1
-
ρ
2
)
-
μ
1
-
γ
1
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
1
d
φ
,
ρ
2
=
ρ
2
(
x
,
0
)
+
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
K
2
(
σ
-
ρ
1
-
ρ
2
)
-
μ
2
-
γ
2
f
Γ
(
α
)
×
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
1
+
k
2
ρ
2
]
d
φ
]
ρ
2
d
φ
.
Since it is not certain that the above has exact solution, the second step of the method proposes that the solution should be in the series of Poincare meaning:
(47)
ρ
1
(
x
,
t
)
=
∑
n
=
0
∞
p
n
ρ
1
n
(
x
,
t
)
,
ρ
2
(
x
,
t
)
=
∑
n
=
0
∞
p
n
ρ
2
n
(
x
,
t
)
.
Then ρ
1,2 can be included in (46) and after including the imbedding parameter p and comparing terms of the same power of p we arrive at the following set of equations:
(48)
ρ
10
(
x
,
t
)
=
ρ
1
(
x
,
0
)
,
ρ
20
(
x
,
t
)
=
ρ
2
(
x
,
0
)
,
ρ
11
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
K
1
(
σ
-
ρ
10
-
ρ
20
)
-
μ
1
-
γ
1
f
Γ
(
α
)
×

∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10

+

k
2
ρ
20
]
d
φ
]
ρ
10
d
ψ
,
ρ
21
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10
+
k
2
ρ
20
]
K
1
(
σ
-
ρ
10
-
ρ
20
)
-
μ
2
-
γ
2
f
Γ
(
α
)
×

∫
0
x
(
x
-
φ
)
α
-
1
[
k
1
ρ
10

+

k
2
ρ
20
]
d
φ
]
ρ
20
d
ψ
.
In general for n greater than 1 we can obtain the rest of the terms with the following iteration formula:
(49)
ρ
1
,
n
+
1
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
K
1
σ
ρ
1
,
n
-
K
1
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
1
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
1
α
n
(
x
,
t
)
]
d
ψ
,
ρ
2
,
n
+
1
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
×
[
K
2
σ
ρ
2
,
n
-
K
2
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
2
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
2
α
n
(
x
,
φ
)
]
d
ψ
.
As soon as one is provided with initial condition, this algorithm can be implemented easily. To make this simple, we provide the summary of the procedure for this case in the following algorithm.Algorithm 2.
(i) Put in:I
1
(
x
)
,
I
2
(
x
)—as preliminary conditions, meaning in our case ρ
1
(
x
,
0
) and ρ
1
(
x
,
0
), k—number terms in the rough calculation.
(ii) Output:ρ
1
approx
(
x
,
t
)
,
ρ
2
approx
(
x
,
t
) the approximate solutions.Step 1.
Putρ
10
(
x
,
t
)
=
I
1
(
x
)
,
ρ
20
(
x
,
t
)
=
I
2
(
x
) and ρ
1
approx
(
x
,
t
)
=
ρ
10
(
x
,
t
)
,
ρ
2
approx
(
x
,
t
)
=
I
2
(
x
).Step 2.
Fork
=
0 to n
-
1 do Step 3, Step 4, and Step 5.Step 3.
Compute(50)
L
1
n
=
K
1
σ
ρ
1
,
n
-
K
1
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
1
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
1
α
n
(
x
,
t
)
,
L
2
n
=
K
2
σ
ρ
2
,
n
-
K
2
∑
j
=
0
n
ρ
1
,
n
-
j
-
1
ρ
1
j
-
K
2
∑
j
=
0
n
ρ
2
,
n
-
j
-
1
ρ
1
j
-
μ
1
ρ
1
,
n
-
ρ
1
n
J
2
α
n
(
x
,
t
)
.Step 4.
Compute(51)
ρ
1
,
n
+
1
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
L
1
n
d
ψ
,
ρ
2
,
n
+
1
(
x
,
t
)
=
∫
0
t
(
t
-
φ
)
α
-
1
Γ
(
α
)
L
2
n
d
ψ
.Step 5.
Computeρ
1
approx
(
x
,
t
)
=
ρ
1
approx
(
x
,
t
)
+
ρ
1
,
n
+
1
(
x
,
t
), ρ
2
approx
(
x
,
t
)
=
ρ
2
approx
(
x
,
t
)
+
ρ
2
,
n
+
1
(
x
,
t
).
Stop.Remark 3.
A direct comparison of both methods shows that they all use integral transform and the Poincare series. However in the first method, one will first apply the Laplace transform and again apply the inverse-Laplace operator which sometime can be very difficult to handle due to the instability of the Laplace transform operator. But the second method uses just a simple integral transform that can be sometime implemented with the new powerful software. We therefore conclude that the second method is much easier to implement rather than the first method. Nevertheless both methods have been proven to be efficient.In the recent decade it was observed by several editors that many scholars using iteration method to find approximate or exact solutions of complicated equations were not investigating the stability and the convergence of their iterations. This observation led us to think that this field was much easier, and other editors said it was a child play. In order to avoid classifying this paper under this class, we will present without loss of generality the convergence of our iteration number 2 and this is done in the next section.


## 4.3. Convergence Investigation

In this section, we examine the convergence property of the approximated solution for the formation of liver zone equation. Let us consider the formation of liver zone equation in the Hilbert spaceH
=
L
2
(
(
η
,
λ
)
×
[
0
,
T
]
) defined as
(52)
H

=


{
(
u
,
v
)
:
(
η
,
λ
)

×

[
0
,
T
]
,
with
∫
(
t
-
φ
)
α
-
1
Γ
(
α
)
u
v
d
ι
d
κ

<

∞
}
.
Then the operator is of the form
(53)
H
(
u
,
v
)
=
u
(
a
(
σ
-
u
-
v
)
-
μ
-
b
∫
0
x
[
c
u
+
d
v
]
d
φ
)
.
That can be reduced to the following operator due to the physical properties of the problem under investigation:
(54)
H
(
u
)
=
u
(
a
(
σ
-
u
-
u
)
-
μ
-
b
∫
0
x
[
c
u
+
d
u
]
d
φ
)
.
The homotopy decomposition method is convergent if the following conditions are satisfiedHypothesis 1.
It is possible for us to find a positive constant, say,F, such that the inner product holds in the Hilbert space H
(55)
(
H
(
u
)
-
H
(
v
)
,
u
-
v
)
≥
F
∥
u
-
v
∥
,
∀
v
,
u
∈
H
.Hypothesis 2.
As far as for allv
,
u
∈
H are bounded this implies, we can find a positive constant, say, C, such that ∥
u
∥
,
∥
v
∥
≤
C, then we can find Φ
(
C
)
>
0 such that
(56)
(
H
(
u
)
-
H
(
v
)
,
g
)
≥
Φ
(
C
)
∥
u
-
v
∥
∥
z
∥
,
∀
z
∈
H
.
We can therefore state the following theorem for the sufficient condition for the convergence of the generalized equation.Theorem 4.
Let us consider(57)
H
(
u
)
=
∂
α
u
∂
t
α
=
u
(
a
(
σ
-
u
-
u
)
-
μ
-
b
∫
0
x
[
c
u
+
d
u
]
d
φ
)
,
w
i
t
h
a
σ
>
2
a
+
b
2
(
c
+
d
)
and consider the initial and boundary condition for the generalized equation; then, the homotopy decomposition method leads to a special solution of system (4).We will present the proof of this theorem by just verifying Hypothesises1 and 2.Proof.
We have that(58)
H
(
u
)
-
H
(
v
)
=
a
σ
(
u
-
v
)
-
2
a
(
u
2
-
v
2
)
-
μ
(
v
-
u
)
-
b
1
u
∫
0
x
[
u
]
d
φ
+
b
1
v
∫
0
x
[
v
]
d
φ
.
The above equation can be roughly converted to
(59)
H
(
u
)
-
H
(
v
)
=
a
σ
(
u
-
v
)
-
2
a
(
u
2
-
v
2
)
-
μ
(
u
-
v
)
-
b
1
2
d
(
u
2
-
v
2
)
d
x
=
a
σ
(
u
-
v
)
-
2
a
(
u
-
v
)
(
u
+
v
)
-
b
1
2
d
[
(
u
-
v
)
(
u
+
v
)
]
d
x
.
However according to the physical property of this problem, we have that ρ
1
+
ρ
2 the total cell density cannot exceed some fixed maximum density σ of cell sites, as limited by familiar phenomenon of contact inhibition; thus, we are urged to conclude that we can find a positive constant, say, h, such that u
+
v
=
h, and then,
(60)
=
a
σ
(
u
-
v
)
-
2
a
(
u
-
v
)
(
u
+
v
)
-
b
1
2
d
[
(
u
-
v
)
(
u
+
v
)
]
d
x
can be converted to
(61)
=
a
σ
(
u
-
v
)
-
2
a
(
u
-
v
)
h
-
b
1
2
d
[
(
u
-
v
)
h
]
d
x
.
With the above in hand, we have the following inner product:
(62)
(
H
(
u
)
-
H
(
v
)
,
u
-
v
)
=
a
σ
(
(
u
-
v
)
,
u
-
v
)
-
2
a
(
(
u
-
v
)
,
u
-
v
)
h
-
h
b
1
2
d
[
(
(
u
-
v
)
,
u
-
v
)
]
d
x
.
According to the physical property of this problem, u
,
v are bounded; therefore, we can find a positive constant F such that (
u
,
u
)
,
(
v
,
v
)
<
F
2. It follows by the use of Schwartz inequality that
(63)
d
[
(
(
u
-
v
)
,
u
-
v
)
]
d
x
≤
∥
(
u
-
v
)
x
∥
∥
u
-
v
∥
.
Now since there exists a positive constant ω such that ∥
(
u
-
v
)
x
∥
≤
ω
∥
u
-
v
∥ then
(64)
d
[
(
(
u
-
v
)
,
u
-
v
)
]
d
x
≤
ω
∥
u
-
v
∥
2
.
On the other hand, we have that
(65)
(
(
u
-
v
)
,
u
-
v
)
=
∥
u
-
v
∥
2
.
Now making a direct substitution of (64) and (65) into (62) we arrive at
(66)
(
H
(
u
)
-
H
(
v
)
,
u
-
v
)
≥
F
∥
u
-
v
∥
2
with F
=
a
σ
-
2
a
h
-
(
h
b
1
/
2
). Therefore Hypothesis 1 is true for the generalized version. We will now present the proof for Hypothesis 2, as we are directly computing
(67)
(
H
(
u
)
-
H
(
v
)
,
z
)
=
a
σ
(
(
u
-
v
)
,
z
)
-

2
a
(
(
u
-
v
)
,
z
)
h

-
h
b
1
2
d
[
(
(
u
-
v
)
,
z
)
]
d
x
≥
Φ
(
C
)
∥
u
-
v
∥
∥
z
∥
.
Therefore we can see that Hypothesis 2 is correct as well. We can conclude without fear that the homotopy decomposition method works perfectly for this generalized equation.


## 5. Conclusion

The fractional order derivative concepts have nowadays gained the world of modelling real world problem due to their properties. This has been confirmed in many studies including almost all branches of sciences. Since the formation of the liver zone is a very sensitive study issue, it has come to our mind to present the mathematical formula underpinning the formation of the liver zone within the scope of fractional derivative. We achieved this by making use of the well-known Caputo derivative because it benefits some physical properties that other including but not limited to the Rieman-Liouville and the Jumarie fractional derivative. We have with care examined the stationary solutions. Since in the recent decade several iteration methods were proposed because formal well-known analytical methods were found to be limited specially when dealing with nonlinear problem, we make use of two different techniques to get to the bottom of the generalized equation. To help readers that are not acquainted with coding theory, we presented the summary of our methods in the form of algorithms. Keeping in mind that one of the greater challenges while using iteration method is to show the convergence, we successfully proved the convergence of the second method for the generalized equation.



---

*Source: 478028-2014-09-08.xml*