# Nonlinear Dynamics and Analysis of Intracranial Saccular Aneurysms with Growth and Remodeling

Authors: Manal Badgaish; Jeng-Eng Lin; Padmanabhan Seshaiyer
Journal: Journal of Nonlinear Dynamics
 (2016)

---

## Abstract

A new mathematical model for the interaction of blood flow with the arterial wall surrounded by cerebral spinal fluid is developed with applications to intracranial saccular aneurysms. The blood pressure acting on the inner arterial wall is modeled via a Fourier series, the arterial wall is modeled as a spring-mass system incorporating growth and remodeling, and the surrounding cerebral spinal fluid is modeled via a simplified one-dimensional compressible Euler equation with inviscid flow and negligible nonlinear effects. The resulting nonlinear coupled fluid-structure interaction problem is analyzed and a perturbation technique is employed to derive the first-order approximation solution to the system. An analytical solution is also derived for the linearized version of the problem using Laplace transforms. The solutions are validated against related work from the literature and the results suggest the biological significance of the inclusion of the growth and remodeling effects on the rupture of intracranial aneurysms.


---

## Body

## 1. Introduction

Intracranial saccular aneurysm is a focal dilatation of the arterial wall that can be found in the circle of Willis which is a network of vessels at the base of the brain. The aneurysm which is a soft tissue interacts with a variety of flows including the blood as well as the cerebral spinal fluid. Based on the influence of various biomechanical factors, the growing aneurysm can potentially rupture that leads to either a neurological disorder or death. About 80% to 90% of ruptured aneurysms lead to death [1]. One of the invasive treatments for intracranial saccular aneurysms is to employ surgical procedures which may be associated with risk and death [2].Over the last two decades, there have been several efforts to investigate the genesis of this disease and to develop a way for prediction of rupture through mathematical modeling [3–6]. Different groups of researchers had identified the elastodynamics of the arterial wall interaction with the blood flow to be the main reason for the rupture of an aneurysm [7–9]. A coupled fluid-structure model to understand the elastodynamics better was later studied more extensively [10–13]. These models introduced mathematical models of increasing complexity for intracranial saccular aneurysms that described the coupled interaction between blood, arterial wall, and cerebral spinal fluid (CSF). In [11], the CSF was modeled using simplified Navier-Stokes equations, whereas the arterial wall structure was modeled using a spring-mass system. A Fourier series was used to model the interaction between blood pressure and inner wall. While the model developed yielded good insight into understanding rupture, this model did not accommodate growth and remodeling effects of the soft-tissue wall.There are three main constituents of the artery wall, namely, the elastin, the collagen, and the smooth muscle [14, 15]. The elastin is a stable protein and is considered the most load bearing element that functions as resistance to the formation of an aneurysm, whereas the collagen is the protein that is responsible for preventing rupture after formation of an aneurysm. The growth of the aneurysm is associated with deficiency of elastin and weakening of the artery wall [16]. Hence, elastin and collagen should be incorporated into the modeling of arterial wall in order to obtain an accurate biological model of the aneurysm that maybe can lead to better interpretation and prediction for this disease. This is one of the main contributions of this work. In the model we will consider in this paper we will assume elastin and collagen to be one-dimensional passive fibers and also, for simplicity, we will not consider the effects of the smooth muscle for this work.In Section2, we will describe the mathematical model that we will consider to solve a coupled fluid-structure problem. Section 3 discusses the nonlinear analysis and perturbation solution for the model. In Section 4, the linearization of the growth and remodeling model is considered and solved analytically using Laplace transformation. After that, some important computational results and comparisons are discussed in Section 5. Finally, we conclude by some biological interpretation of the results and suggestion for future work.


## 2. Models and Background

The mathematical modeling of intracranial saccular aneurysm builds on models developed in [11, 12] which included a simple mathematical model of a thin-walled, spherical intracranial aneurysm surrounded by cerebral spinal fluid which is referred to as CSF (see Figure 1).Figure 1
Model of an aneurysm in an arterial wall with blood inside and CSF outside.The model was based on three separate components: modeling the blood pressure acting on the inner wall of the arteries, modeling the CSF that is located outside the wall, and modeling the arterial wall itself. Combining the various components yields a coupled fluid-structure interaction model. This model assumed a simple spring-mass model for the wall and the focus of this paper is to enhance this wall model to incorporate realistic nonlinear soft-tissue behavior.


### 2.1. Model of the Blood Pressure

The blood pressure is modeled using finite Fourier series since we consider the behavior to be pulsatile [10, 17, 18]:(1)
P
B
t
=
P
m
+
∑
n
=
1
N
A
n
cos
⁡
n
w
t
+
B
n
sin
⁡
n
w
t
,where P
m is the mean blood pressure, w is the fundamental circular frequency, and A
n and B
n are the Fourier coefficients for N harmonic.


### 2.2. Model of the Arterial Wall

We consider the arterial wall to be modeled using a simple spring-mass system that incorporates the elastin and collagen effects in the outer wall of the arteries. The force of this system may be denoted byF
S which is given by(2)
F
S
=
K
E
A
E
t
σ
E
ϵ
E
+
K
C
A
C
t
σ
C
ϵ
C
-
a
P
B
t
,where K
E and  K
C are the scaling coefficients, A
E
(
t
) and A
C
(
t
) are the cross-sectional areas, and σ
E
(
ϵ
E
) and σ
C
(
ϵ
C
) are the stresses for elastin and collagen, respectively. F
S may be considered the force exerted on the outer wall. It is expressed in terms of the respective forces in the elastin and collagen along with the force from the blood pressure in the opposite direction. Each of the respective forces from the elastin and collagen is proportional to the respective stresses σ
E and σ
C. These forces help to balance the force expressed by the blood flow on the inner wall.In this work we assume that these stresses are related to the respective strains through linear constitutive laws (material linearity) given by(3)
σ
E
=
ϵ
E
,
σ
C
=
ϵ
Cand the strains are assumed to be nonlinear (geometric nonlinearity) given by(4)
ϵ
E
=
2
L
u
+
u
2
2
L
2
,
ϵ
C
=
2
L
u
+
u
2
+
L
2
1
-
r
2
2
L
2
r
2
,where L denotes the length of the unstrained tissue, u its extension, and r is the stretched factor of unstrained tissue of collagen fiber [6].


### 2.3. Model of the Cerebral Spinal Fluid

The model of CSF that is considered in this paper is the simplified one-dimensional compressible Euler equation with inviscid flow. Assuming negligible nonlinear effects and a constitutive relationship for the pressure [11], one can derive the following wave equation:(5)
v
t
=
c
2
u
x
x
,
u
t
=
v
.Here u
(
x
,
t
) is assumed to be the displacement of the CSF with v
(
x
,
t
) as the velocity. Since we are looking to find the movement of outer wall due to the interaction with CSF, we consider x
=
0 to denote the outer wall and therefore we are interested in finding the solution to (5) at x
=
0 that will describe the movement of the wall at any time t
≥
0. In order to solve the system, we will assume that the displacement and velocity of the CSF is zero initially. This is given by the initial conditions:(6)
u
x
,
0
=
v
x
,
0
=
0
.The boundary conditions will be described later after the discussion of the modeling of the blood pressure and the arterial wall.


### 2.4. Governing Equation of Motion

In order to solve system (5), we need two boundary conditions. The first boundary condition is at point x
=
0, and it can be derived from the model of blood pressure and the arterial wall that we have discussed. Note that the force balance equation at x
=
0 may be written as(7)
F
T
=
F
F
-
F
S
.Here,(8)
F
T
=
m
u
t
t
0
,
tis the total force, where m is the mass of the wall. Consider that(9)
F
F
=
-
P
a
=
ρ
c
2
u
x
0
,
t
ais the force of the fluid which is the product of the pressure P and the cross-sectional area a. The pressure is given by P
=
-
ρ
c
2
u
x, where ρ is the density of the CSF which is assumed to be slightly compressible [11]. The total force of the spring from (2) is given by(10)
F
S
=
K
E
A
E
ϵ
E
+
K
C
A
C
ϵ
C
-
a
P
B
t
.Then, substituting (8), (9), and (10) into (7) we obtain the first boundary condition at x
=
0:(11)
m
u
t
t
0
,
t
=
ρ
c
2
a
u
x
0
,
t
-
K
E
A
E
ϵ
E
-
K
C
A
C
ϵ
C
+
a
P
m
+
∑
n
=
1
N
a
A
n
cos
⁡
n
w
t
+
a
B
n
sin
⁡
n
w
t
.The second boundary condition can be obtained using the plane wave approximation that states that the waves from the wall will die down some fixed distance away from the wall. If this can be applied at point x
=
L, then the second boundary condition becomes [11](12)
v
L
,
t
=
-
c
u
x
L
,
t
.Combining (5), (6), (11), and (12), we obtain the following coupled fluid-structure interaction problem:(13)
u
t
t
=
c
2
u
x
x
,
u
x
,
0
=
u
t
x
,
0
=
0
,
m
u
t
t
0
,
t
=
a
P
B
t
+
ρ
c
2
a
u
x
0
,
t
+
-
K
E
A
E
L
-
K
C
A
C
L
r
2
u
0
,
t
+
-
K
E
A
E
2
L
2
-
K
C
A
C
2
L
2
r
2
u
2
0
,
t
-
K
C
A
C
1
-
r
2
2
r
2
,
u
t
L
,
t
=
-
c
u
x
L
,
t
.Solving the above coupled system will yield u
(
x
,
t
) which is the displacement of the fluid in the domain [
0
,
L
]. Note that the displacement of the outer boundary of the wall is the same as u
(
0
,
t
) which is the displacement of the fluid at the left end boundary of its domain. As described in [11], we will also study the movement of the outer wall, namely, u
(
0
,
t
), to be an indicator of the potential for rupture of the aneurysm. The higher the value of u
(
0
,
t
), the more the potential for rupture. However, one must solve nonlinear coupled system (13) which does not admit an exact solution. One must therefore apply other approaches to obtain this solution. Next, we will employ a nonlinear analysis and perturbation method in this model to obtain an approximate solution to system (13) which is one of the central contributions of this work.


## 2.1. Model of the Blood Pressure

The blood pressure is modeled using finite Fourier series since we consider the behavior to be pulsatile [10, 17, 18]:(1)
P
B
t
=
P
m
+
∑
n
=
1
N
A
n
cos
⁡
n
w
t
+
B
n
sin
⁡
n
w
t
,where P
m is the mean blood pressure, w is the fundamental circular frequency, and A
n and B
n are the Fourier coefficients for N harmonic.


## 2.2. Model of the Arterial Wall

We consider the arterial wall to be modeled using a simple spring-mass system that incorporates the elastin and collagen effects in the outer wall of the arteries. The force of this system may be denoted byF
S which is given by(2)
F
S
=
K
E
A
E
t
σ
E
ϵ
E
+
K
C
A
C
t
σ
C
ϵ
C
-
a
P
B
t
,where K
E and  K
C are the scaling coefficients, A
E
(
t
) and A
C
(
t
) are the cross-sectional areas, and σ
E
(
ϵ
E
) and σ
C
(
ϵ
C
) are the stresses for elastin and collagen, respectively. F
S may be considered the force exerted on the outer wall. It is expressed in terms of the respective forces in the elastin and collagen along with the force from the blood pressure in the opposite direction. Each of the respective forces from the elastin and collagen is proportional to the respective stresses σ
E and σ
C. These forces help to balance the force expressed by the blood flow on the inner wall.In this work we assume that these stresses are related to the respective strains through linear constitutive laws (material linearity) given by(3)
σ
E
=
ϵ
E
,
σ
C
=
ϵ
Cand the strains are assumed to be nonlinear (geometric nonlinearity) given by(4)
ϵ
E
=
2
L
u
+
u
2
2
L
2
,
ϵ
C
=
2
L
u
+
u
2
+
L
2
1
-
r
2
2
L
2
r
2
,where L denotes the length of the unstrained tissue, u its extension, and r is the stretched factor of unstrained tissue of collagen fiber [6].


## 2.3. Model of the Cerebral Spinal Fluid

The model of CSF that is considered in this paper is the simplified one-dimensional compressible Euler equation with inviscid flow. Assuming negligible nonlinear effects and a constitutive relationship for the pressure [11], one can derive the following wave equation:(5)
v
t
=
c
2
u
x
x
,
u
t
=
v
.Here u
(
x
,
t
) is assumed to be the displacement of the CSF with v
(
x
,
t
) as the velocity. Since we are looking to find the movement of outer wall due to the interaction with CSF, we consider x
=
0 to denote the outer wall and therefore we are interested in finding the solution to (5) at x
=
0 that will describe the movement of the wall at any time t
≥
0. In order to solve the system, we will assume that the displacement and velocity of the CSF is zero initially. This is given by the initial conditions:(6)
u
x
,
0
=
v
x
,
0
=
0
.The boundary conditions will be described later after the discussion of the modeling of the blood pressure and the arterial wall.


## 2.4. Governing Equation of Motion

In order to solve system (5), we need two boundary conditions. The first boundary condition is at point x
=
0, and it can be derived from the model of blood pressure and the arterial wall that we have discussed. Note that the force balance equation at x
=
0 may be written as(7)
F
T
=
F
F
-
F
S
.Here,(8)
F
T
=
m
u
t
t
0
,
tis the total force, where m is the mass of the wall. Consider that(9)
F
F
=
-
P
a
=
ρ
c
2
u
x
0
,
t
ais the force of the fluid which is the product of the pressure P and the cross-sectional area a. The pressure is given by P
=
-
ρ
c
2
u
x, where ρ is the density of the CSF which is assumed to be slightly compressible [11]. The total force of the spring from (2) is given by(10)
F
S
=
K
E
A
E
ϵ
E
+
K
C
A
C
ϵ
C
-
a
P
B
t
.Then, substituting (8), (9), and (10) into (7) we obtain the first boundary condition at x
=
0:(11)
m
u
t
t
0
,
t
=
ρ
c
2
a
u
x
0
,
t
-
K
E
A
E
ϵ
E
-
K
C
A
C
ϵ
C
+
a
P
m
+
∑
n
=
1
N
a
A
n
cos
⁡
n
w
t
+
a
B
n
sin
⁡
n
w
t
.The second boundary condition can be obtained using the plane wave approximation that states that the waves from the wall will die down some fixed distance away from the wall. If this can be applied at point x
=
L, then the second boundary condition becomes [11](12)
v
L
,
t
=
-
c
u
x
L
,
t
.Combining (5), (6), (11), and (12), we obtain the following coupled fluid-structure interaction problem:(13)
u
t
t
=
c
2
u
x
x
,
u
x
,
0
=
u
t
x
,
0
=
0
,
m
u
t
t
0
,
t
=
a
P
B
t
+
ρ
c
2
a
u
x
0
,
t
+
-
K
E
A
E
L
-
K
C
A
C
L
r
2
u
0
,
t
+
-
K
E
A
E
2
L
2
-
K
C
A
C
2
L
2
r
2
u
2
0
,
t
-
K
C
A
C
1
-
r
2
2
r
2
,
u
t
L
,
t
=
-
c
u
x
L
,
t
.Solving the above coupled system will yield u
(
x
,
t
) which is the displacement of the fluid in the domain [
0
,
L
]. Note that the displacement of the outer boundary of the wall is the same as u
(
0
,
t
) which is the displacement of the fluid at the left end boundary of its domain. As described in [11], we will also study the movement of the outer wall, namely, u
(
0
,
t
), to be an indicator of the potential for rupture of the aneurysm. The higher the value of u
(
0
,
t
), the more the potential for rupture. However, one must solve nonlinear coupled system (13) which does not admit an exact solution. One must therefore apply other approaches to obtain this solution. Next, we will employ a nonlinear analysis and perturbation method in this model to obtain an approximate solution to system (13) which is one of the central contributions of this work.


## 3. Nonlinear Dynamics and Analysis

Nonlinear analysis and perturbation method techniques are powerful tools to study the behavior of complex nonlinear models and help to obtain an approximate solution to these models. In this work, the nonlinear phase analysis and perturbation technique for (13) will be presented. We will employ the following steps. First we will convert (13) which is an initial-boundary value partial differential equation system into an initial value ordinary differential equation system. Following this we will perform a phase analysis of the newly obtained initial value differential equation. We will then apply the perturbation techniques to the resulting equations. These steps are described next.


### 3.1. Converting the PDE System into Initial Value Differential Equation

Since the governing equation of system (13) is a one-dimensional wave equation, it has a unique solution of the form [19](14)
u
x
,
t
=
f
~
x
-
c
t
+
g
~
x
+
c
t
.Differentiating the above solution with respect to variables x and t yields(15)
u
t
x
,
t
=
-
c
f
~
′
x
-
c
t
+
c
g
~
′
x
+
c
t
,
(16)
u
t
t
x
,
t
=
c
2
f
~
′
′
x
-
c
t
+
c
2
g
~
′
′
x
+
c
t
,
(17)
u
x
x
,
t
=
f
~
′
x
-
c
t
+
g
~
′
x
+
c
t
,
(18)
u
x
x
x
,
t
=
f
~
′
′
x
-
c
t
+
g
~
′
′
x
+
c
t
.Substituting (15)–(17) in the boundary condition at x
=
L of system (13) yields(19)
-
c
f
~
′
L
-
c
t
+
c
g
~
′
L
+
c
t
=
-
c
f
~
′
x
-
c
t
-
c
g
~
′
x
+
c
t
.This yields g
~
′
(
L
+
c
t
)
=
0 and therefore(20)
g
~
L
+
c
t
=
D
,where D is arbitrary constant. Using (14), the unique solution of the wave equation becomes (21)
u
x
,
t
=
f
~
x
-
c
t
+
D
.Since D is an arbitrary constant, it can be absorbed inside the function f
~
(
x
-
c
t
). Therefore, the general solution of the system is(22)
u
x
,
t
=
f
x
-
c
t
.Since our goal is to determine the solution u
(
0
,
t
), next we substitute (22) into the boundary condition at x
=
0 of system (13) which yields(23)
m
c
2
f
′
′
-
c
t
=
a
P
B
t
+
ρ
c
2
a
f
′
-
c
t
+
-
K
E
A
E
L
-
K
C
A
C
L
r
2
f
-
c
t
+
-
K
E
A
E
2
L
2
-
K
C
A
C
2
L
2
r
2
f
2
-
c
t
-
K
C
A
C
1
-
r
2
2
r
2which can be rewritten as(24)
m
c
2
f
′
′
-
c
t
-
ρ
c
2
a
f
′
-
c
t
+
D
1
f
-
c
t
+
D
2
f
2
-
c
t
=
a
P
B
t
+
D
3with the following initial condition: f
(
0
)
=
0, f
′
(
0
)
=
0. Here, (25)
P
B
t
=
P
m
+
∑
n
=
1
N
A
n
cos
⁡
n
w
t
+
B
n
sin
⁡
n
w
t
,
D
1
=
K
E
A
E
L
+
K
C
A
C
L
r
2
,
D
2
=
K
E
A
E
2
L
2
+
K
C
A
C
2
L
2
r
2
,
D
3
=
-
K
C
A
C
1
-
r
2
2
r
2
.Lemma 1.
Nonlinear differential equation (24) can be rewritten as the following nonlinear second-order differential equation:(26)
m
w
2
h
′
′
t
+
ρ
c
w
a
h
′
t
+
D
1
h
t
+
D
2
h
2
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
cos
⁡
n
t
+
B
n
sin
⁡
n
t
.Proof.
To prove this, we will use a change of variable. Let us defineT
=
-
c
t and substitute this in differential equation (24) to obtain
(27)
m
c
2
f
′
′
T
-
ρ
c
2
a
f
′
T
+
D
1
f
T
+
D
2
f
2
T
=
a
P
m
+
a
∑
n
=
1
N
A
n
cos
⁡
n
w
T
-
c
+
B
n
sin
⁡
n
w
T
-
c
+
D
3
.Next, we will choose another change of variable and let t
~
=
w
T
/
-

c. Also define a new function h
(
t
~
)
=
f
(
T
). Then applying chain rule, we have(28)
h
′
t
~
=
-
c
w
f
′
T
,
h
′
′
t
~
=
-
c
w
2
f
′
′
T
.Using these relations and renaming t
~
=
t, (27) can be rewritten as(29)
m
w
2
h
′
′
t
+
ρ
c
w
a
h
′
t
+
D
1
h
t
+
D
2
h
2
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
cos
⁡
n
t
+
B
n
sin
⁡
n
t
.Note that the last summation term in (29) is now 2
π periodic.


### 3.2. Phase Analysis of the Initial Value Differential Equation

Consider the homogeneous case for (29). Dividing throughout by m
w
2 yields(30)
h
′
′
t
+
Q
1
h
′
t
+
Q
2
h
t
+
Q
3
h
2
t
=
0
,where(31)
Q
1
=
ρ
c
a
m
w
,
Q
2
=
D
1
m
w
2
,
Q
3
=
D
2
m
w
2
.Defining a new vector,(32)
H
t
=
h
~
1
t
h
~
2
t
,where h
~
1
(
t
)
=
h
(
t
) and h
~
2
(
t
)
=
h
′
(
t
); we can now rewrite (30) into two first-order equations as follows:(33)
h
~
1
′
=
h
~
2
,
h
~
2
′
=
-
Q
1
h
~
2
-
Q
2
h
~
1
-
Q
3
h
~
1
2
.It is trivial to note from the system above that there are two equilibrium points, namely, the origin (
0,0
) and the point -
Q
2
/
Q
3
,
0.Considering the first equilibrium point, that is(
0,0
), system (33) can be rewritten in a matrix form as(34)
H
′
t
=
A
H
t
+
g
t
,where(35)
H
′
t
=
h
~
1
′
h
~
2
′
,
A
=
D
h
~
0
,
0
=
0
1
-
Q
2
-
Q
1
,
h
~
=
h
~
2
-
Q
1
h
~
2
-
Q
2
h
~
1
-
Q
3
h
~
1
2
,
H
t
=
h
~
1
h
~
2
,
g
t
=
0
-
Q
3
h
~
1
2
t
.The linearized part of (34) is(36)
H
′
t
=
A
H
t
.The eigenvalues of the matrix A are two negative real numbers that can be calculated using(37)
λ
1,2
=
-
Q
1
±
Q
1
2
-
4
Q
2
2and since none of the eigenvalues of A has zero real part, the equilibrium point is hyperbolic. Hence, the behavior of nonlinear system (34) is topologically equivalent to the behavior of linearized system (36), and the approximate solution for this linear system will be good first approximation for nonlinear system (34) [20]. Note that we will restrict the values of ρ,  c,  a,  m,  w,  K
E,  K
C,  A
E,  A
C,  L, and r in order to have two distinct real eigenvalues. This means that(38)
Q
1
2
>
4
Q
2
.This condition as we will see later has biological significance in relation to the displacement of the outer wall. Analyzing the behavior of linear system, it can be seen that the origin (
0,0
) is a stable node (sink) for linear system (36). Thus, it is a stable node for nonlinear system (34) as well [20]. Therefore, in area around the origin the solution is supposed to approach the equilibrium point (see Figure 2).Figure 2
Phase portrait of the system.In order to study the behavior of the nonlinear solution near the second equilibrium point, that is(
-
Q
2
/
Q
3
,
0
), this equilibrium point has to be translated to the origin using the following transformation:(39)
h
1
^
=
h
~
1
+
Q
2
Q
3
,
h
2
^
=
h
~
2
.Using above transformation, the new system can be rewritten as(40)
h
1
^
′
=
h
2
^
,
h
2
^
′
=
-
Q
1
h
2
^
+
Q
2
h
1
^
-
Q
3
h
1
^
2
.This can be written in matrix form as(41)
H
t
^
′
=
A
^
H
t
^
+
g
t
^
,where(42)
H
t
^
′
=
h
1
^
′
h
2
^
′
,
A
^
=
D
h
^
0
,
0
=
0
1
Q
2
-
Q
1
,
h
^
=
h
2
^
-
Q
1
h
2
^
+
Q
2
h
1
^
-
Q
3
h
1
^
2
,
H
t
^
=
h
1
^
h
2
^
,
g
t
^
=
0
-
Q
3
h
1
^
2
t
.The linearized part of (41) is(43)
H
t
^
′
=
A
^
H
t
^
.The eigenvalues of the matrix A
^ are two distinct real numbers with opposite signs and can be calculated using(44)
λ
1,2
=
-
Q
1
±
Q
1
2
+
4
Q
2
2
.Therefore, the behavior of nonlinear system (41) is topologically equivalent to the behavior of linearized system (43). We have found that this equilibrium point -
Q
2
/
Q
3
,
0 is unstable saddle for the linear system; therefore, it is also an unstable saddle for the nonlinear system [20].


### 3.3. Existence of the Solution to Perturbed Nonlinear Nonhomogeneous System

Let us consider again the following second-order nonhomogeneous weakly nonlinear ODE:(45)
h
′
′
t
+
Q
1
h
′
t
+
Q
2
h
t
+
ϵ
Q
3
h
2
t
=
F
t
,where(46)
F
t
=
Q
4
a
P
m
+
D
3
+
a
Q
4
∑
n
=
1
N
A
n
cos
⁡
n
t
+
B
n
sin
⁡
n
tis 2
π periodic Fourier series and Q
4
=
1
/
m
w
2. This system can be written in vector form as(47)
h
~
1
′
=
h
~
2
,
h
~
2
′
=
-
Q
1
h
~
2
-
Q
2
h
~
1
+
F
t
-
ϵ
Q
3
h
~
1
2
,where h
~
1
(
t
)
=
h
(
t
) and h
~
2
(
t
)
=
h
′
(
t
).The associated reduced system whenϵ
=
0 is(48)
h
~
1
′
=
h
~
2
,
h
~
2
′
=
-
Q
1
h
~
2
-
Q
2
h
~
1
+
F
t
.A particular solution of reduced system (48) is expressed as(49)
ψ
t
=
Q
4
a
P
m
+
D
3
Q
2
0
+
∑
n
=
1
N
1
Q
2
-
n
2
2
+
n
Q
1
2
a
n
∗
n
b
n
∗
cos
⁡
n
t
+
b
n
∗
-
n
a
n
∗
sin
⁡
n
t
,where(50)
a
n
∗
=
Q
2
-
n
2
Q
4
A
n
-
n
Q
1
Q
4
B
n
,
b
n
∗
=
n
Q
1
Q
4
A
n
+
Q
2
-
n
2
Q
4
B
n
.Note that, from here, we can obtain all the general solutions for (48) which will require the application of the initial conditions h
~
1
(
0
)
=
h
~
2
(
0
)
=
0 to solve for any unknown constants in the general solution. Then, (45) has a 2
π periodic solution h
=
h
(
t
,
ϵ
) that is analytic in ϵ. Therefore, the solution can be expressed as convergent power series [21]:(51)
h
t
,
ϵ
=
h
0
t
+
ϵ
h
1
t
+
ϵ
2
h
2
t
+
⋯
=
∑
n
=
0
∞
ϵ
n
h
n
t
,where each function h
n
(
t
) is 2
π periodic in t.


### 3.4. Perturbation Technique to Find First-Order Approximation

Let us consider (24) with perturbed term:(52)
m
c
2
f
′
′
-
c
t
-
ρ
c
2
a
f
′
-
c
t
+
D
1
f
-
c
t
+
ϵ
D
2
f
2
-
c
t
=
F
twith the following initial condition: f
(
0
)
=
0, f
′
(
0
)
=
0, where(53)
F
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
cos
⁡
n
w
t
+
B
n
sin
⁡
n
w
tand ϵ here is assumed to be very small. Following the ideas in Lemma 1, the above differential equation can be rewritten as(54)
m
w
2
h
′
′
t
+
ρ
c
w
a
h
′
t
+
D
1
h
t
+
ϵ
D
2
h
2
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
cos
⁡
n
t
+
B
n
sin
⁡
n
twith h
(
0
)
=
0
,
h
′
(
0
)
=
0. Let the solution of differential equation (54) be expressed as power series as in the previous section. Replacing the function h
(
t
) by its power series (51) and comparing coefficients of each order, we obtain the following [22–24].Zeroth-Order Approximation. Consider(55)
m
w
2
h
0
′
′
t
+
ρ
c
w
a
h
0
′
t
+
D
1
h
0
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
c
o
s
n
t
+
B
n
s
i
n
n
t
,
h
0
0
=
0
,
h
0
′
0
=
0
.First-Order Approximation. Consider(56)
m
w
2
h
1
′
′
t
+
ρ
c
w
a
h
1
′
t
+
D
1
h
1
t
+
D
2
h
0
t
2
=
0
,
h
1
0
=
0
,
h
1
′
0
=
0
.We can now solve (55) that yields(57)
h
0
t
=
B
~
1
e
λ
1
t
+
B
~
2
e
λ
2
t
+
a
0
+
∑
n
=
1
N
a
n
cos
⁡
n
t
+
b
n
sin
⁡
n
t
,where (58)
λ
1,2
=
-
ρ
c
w
a
±
ρ
c
w
a
2
-
4
m
w
2
D
1
2
m
w
2
,
a
0
=
a
P
m
+
D
3
D
1
,
a
n
=
a
D
1
-
n
2
w
2
m
A
n
-
a
n
ρ
w
a
c
B
n
D
1
-
n
2
w
2
m
2
+
n
ρ
w
a
c
2
,
b
n
=
a
n
ρ
w
a
c
A
n
+
a
D
1
-
n
2
w
2
m
B
n
D
1
-
n
2
w
2
m
2
+
n
ρ
w
a
c
2
,
B
~
1
=
λ
2
a
0
+
∑
n
=
0
N
a
n
-
∑
n
=
1
N
n
b
n
λ
1
-
λ
2
,
B
~
2
=
λ
1
a
0
+
∑
n
=
0
N
a
n
-
∑
n
=
1
N
n
b
n
λ
2
-
λ
1
.Since we have restriction condition from Section 3.2, that is, Q
1
2
>
4
Q
2, this is identical to the following condition:(59)
ρ
c
w
a
2
>
4
m
w
2
D
1
.Hence, λ
1
≠
λ
2
.Using solutionh
0
(
t
) from (57), we can rewrite (56) in the form(60)
m
w
2
h
1
′
′
t
+
ρ
c
w
a
h
1
′
t
+
D
1
h
1
t
+
2
D
2
B
~
1
∑
n
=
0
N
a
n
cos
⁡
n
t
e
λ
1
t
+
2
D
2
B
~
1
∑
n
=
1
N
b
n
sin
⁡
n
t
e
λ
1
t
+
2
D
2
B
~
2
∑
n
=
0
N
a
n
cos
⁡
n
t
e
λ
2
t
+
2
D
2
B
~
2
∑
n
=
1
N
b
n
sin
⁡
n
t
e
λ
2
t
+
D
2
B
~
1
2
e
2
λ
1
t
+
2
D
2
B
~
1
B
~
2
e
λ
1
+
λ
2
t
+
D
2
B
~
2
2
e
2
λ
2
t
+
D
2
a
0
2
+
2
D
2
a
0
∑
n
=
1
N
a
n
cos
⁡
n
t
+
b
n
sin
⁡
n
t
+
D
2
2
∑
n
=
1
N
a
n
2
+
D
2
2
∑
n
=
1
N
b
n
2
+
D
2
2
∑
n
=
0
N
a
n
2
cos
⁡
2
n
t
+
D
2
∑
n
=
1
N
a
n
b
n
sin
⁡
2
n
t
-
D
2
2
∑
n
=
0
N
b
n
2
cos
⁡
2
n
t
=
0
.Solution to (60) can be derived to be(61)
h
1
t
=
X
1
+
S
0
e
λ
1
t
+
X
2
+
W
0
e
λ
2
t
+
U
0
+
E
0
+
∑
n
=
1
N
U
n
cos
⁡
n
t
+
∑
n
=
1
N
Q
n
sin
⁡
n
t
+
Z
1
e
2
λ
1
t
+
Z
2
e
2
λ
2
t
+
Z
3
e
λ
1
+
λ
2
t
+
∑
n
=
1
N
E
n
cos
⁡
2
n
t
+
∑
n
=
1
N
G
n
sin
⁡
2
n
t
+
∑
n
=
1
N
S
n
cos
⁡
n
t
e
λ
1
t
+
∑
n
=
1
N
V
n
sin
⁡
n
t
e
λ
1
t
+
∑
n
=
1
N
W
n
cos
⁡
n
t
e
λ
2
t
+
∑
n
=
1
N
Y
n
sin
⁡
n
t
e
λ
2
t
,where (62)
U
0
=
-
D
2
a
0
+
∑
n
=
1
N
a
n
2
/
2
+
∑
n
=
1
N
b
n
2
/
2
D
1
,
S
0
=
-
2
D
2
B
~
1
a
0
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
,
E
0
=
-
D
2
a
0
2
2
D
1
,
W
0
=
-
2
D
2
B
~
2
a
0
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
,
U
n
=
-
2
D
2
a
0
D
1
-
n
2
m
w
2
a
n
+
ρ
c
a
n
w
b
n
D
1
-
n
2
m
w
2
2
+
ρ
c
a
n
w
2
,
Q
n
=
2
D
2
a
0
ρ
c
a
n
w
a
n
-
D
1
-
n
2
m
w
2
b
n
D
1
-
n
2
m
w
2
2
+
ρ
c
a
n
w
2
,
Z
1
=
-
D
2
B
~
1
2
4
m
w
2
λ
1
2
+
2
ρ
c
w
a
λ
1
+
D
1
,
Z
2
=
-
D
2
B
~
2
2
4
m
w
2
λ
2
2
+
2
ρ
c
w
a
λ
1
+
D
1
,
Z
3
=
-
2
D
2
B
~
1
B
~
2
m
w
2
λ
1
+
λ
2
2
+
ρ
c
w
a
λ
1
+
λ
2
+
D
1
,
E
n
=
D
2
D
1
-
4
m
n
2
w
2
/
2
b
n
2
-
a
n
2
+
2
ρ
c
a
n
w
a
n
b
n
D
1
-
4
m
n
2
w
2
2
+
2
ρ
c
a
n
w
2
,
G
n
=
-
D
2
2
ρ
c
a
n
w
/
2
a
n
2
-
b
n
2
+
D
1
-
4
m
n
2
w
2
a
n
b
n
D
1
-
4
m
n
2
w
2
2
+
2
ρ
c
a
n
w
2
,
S
n
=
2
D
2
B
~
1
2
n
w
2
m
λ
1
+
n
ρ
w
a
c
b
n
-
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
-
m
n
2
w
2
a
n
2
n
w
2
m
λ
1
+
ρ
c
a
n
w
2
+
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
-
m
n
2
w
2
2
,
V
n
=
-
2
D
2
B
~
1
2
n
w
2
m
λ
1
+
n
ρ
w
c
a
a
n
+
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
-
m
n
2
w
2
b
n
2
n
w
2
m
λ
1
+
ρ
c
a
n
w
2
+
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
-
m
n
2
w
2
2
,
W
n
=
2
D
2
B
~
2
2
n
w
2
m
λ
2
+
ρ
n
w
m
c
b
n
-
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
-
m
n
2
w
2
a
n
2
n
w
2
m
λ
2
+
ρ
c
a
n
w
2
+
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
-
m
n
2
w
2
2
,
Y
n
=
-
2
D
2
B
~
2
2
n
w
m
λ
2
+
n
ρ
w
m
c
a
n
+
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
-
m
n
2
w
2
b
n
2
n
w
2
m
λ
2
+
ρ
c
a
n
w
2
+
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
-
m
n
2
w
2
2
,
X
1
=
P
1
+
Φ
1
λ
1
-
λ
2
,
X
2
=
P
2
+
Φ
2
λ
2
-
λ
1
,where (63)
P
1
=
λ
2
-
2
λ
1
Z
1
-
λ
2
Z
2
-
λ
1
Z
3
+
λ
2
U
0
,
Φ
1
=
λ
2
∑
n
=
1
N
U
n
-
∑
n
=
1
N
2
n
G
n
+
λ
2
-
λ
1
∑
n
=
0
N
S
n
+
λ
2
∑
n
=
0
N
E
n
-
∑
n
=
1
N
n
V
n
+
Y
n
+
Q
n
,
P
2
=
-
λ
1
Z
1
+
λ
1
-
2
λ
2
Z
2
-
λ
2
Z
3
+
λ
1
U
0
,
Φ
2
=
λ
1
∑
n
=
1
N
U
n
-
∑
n
=
1
N
2
n
G
n
+
λ
1
-
λ
2
∑
n
=
0
N
W
n
+
λ
1
∑
n
=
0
N
E
n
-
∑
n
=
1
N
n
V
n
+
Y
n
+
Q
n
.For simplicity, we restricted the solution to be first-order approximation, so(64)
h
t
=
h
0
t
+
ϵ
h
1
tand then we can have a first-order approximation of the solution we are looking for which is given by(65)
u
0
,
t
=
f
-
c
t
=
f
T
=
h
t
=
h
0
t
+
ϵ
h
1
t
,where h
0 and  h
1 are defined in (57) and (61), respectively.


## 3.1. Converting the PDE System into Initial Value Differential Equation

Since the governing equation of system (13) is a one-dimensional wave equation, it has a unique solution of the form [19](14)
u
x
,
t
=
f
~
x
-
c
t
+
g
~
x
+
c
t
.Differentiating the above solution with respect to variables x and t yields(15)
u
t
x
,
t
=
-
c
f
~
′
x
-
c
t
+
c
g
~
′
x
+
c
t
,
(16)
u
t
t
x
,
t
=
c
2
f
~
′
′
x
-
c
t
+
c
2
g
~
′
′
x
+
c
t
,
(17)
u
x
x
,
t
=
f
~
′
x
-
c
t
+
g
~
′
x
+
c
t
,
(18)
u
x
x
x
,
t
=
f
~
′
′
x
-
c
t
+
g
~
′
′
x
+
c
t
.Substituting (15)–(17) in the boundary condition at x
=
L of system (13) yields(19)
-
c
f
~
′
L
-
c
t
+
c
g
~
′
L
+
c
t
=
-
c
f
~
′
x
-
c
t
-
c
g
~
′
x
+
c
t
.This yields g
~
′
(
L
+
c
t
)
=
0 and therefore(20)
g
~
L
+
c
t
=
D
,where D is arbitrary constant. Using (14), the unique solution of the wave equation becomes (21)
u
x
,
t
=
f
~
x
-
c
t
+
D
.Since D is an arbitrary constant, it can be absorbed inside the function f
~
(
x
-
c
t
). Therefore, the general solution of the system is(22)
u
x
,
t
=
f
x
-
c
t
.Since our goal is to determine the solution u
(
0
,
t
), next we substitute (22) into the boundary condition at x
=
0 of system (13) which yields(23)
m
c
2
f
′
′
-
c
t
=
a
P
B
t
+
ρ
c
2
a
f
′
-
c
t
+
-
K
E
A
E
L
-
K
C
A
C
L
r
2
f
-
c
t
+
-
K
E
A
E
2
L
2
-
K
C
A
C
2
L
2
r
2
f
2
-
c
t
-
K
C
A
C
1
-
r
2
2
r
2which can be rewritten as(24)
m
c
2
f
′
′
-
c
t
-
ρ
c
2
a
f
′
-
c
t
+
D
1
f
-
c
t
+
D
2
f
2
-
c
t
=
a
P
B
t
+
D
3with the following initial condition: f
(
0
)
=
0, f
′
(
0
)
=
0. Here, (25)
P
B
t
=
P
m
+
∑
n
=
1
N
A
n
cos
⁡
n
w
t
+
B
n
sin
⁡
n
w
t
,
D
1
=
K
E
A
E
L
+
K
C
A
C
L
r
2
,
D
2
=
K
E
A
E
2
L
2
+
K
C
A
C
2
L
2
r
2
,
D
3
=
-
K
C
A
C
1
-
r
2
2
r
2
.Lemma 1.
Nonlinear differential equation (24) can be rewritten as the following nonlinear second-order differential equation:(26)
m
w
2
h
′
′
t
+
ρ
c
w
a
h
′
t
+
D
1
h
t
+
D
2
h
2
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
cos
⁡
n
t
+
B
n
sin
⁡
n
t
.Proof.
To prove this, we will use a change of variable. Let us defineT
=
-
c
t and substitute this in differential equation (24) to obtain
(27)
m
c
2
f
′
′
T
-
ρ
c
2
a
f
′
T
+
D
1
f
T
+
D
2
f
2
T
=
a
P
m
+
a
∑
n
=
1
N
A
n
cos
⁡
n
w
T
-
c
+
B
n
sin
⁡
n
w
T
-
c
+
D
3
.Next, we will choose another change of variable and let t
~
=
w
T
/
-

c. Also define a new function h
(
t
~
)
=
f
(
T
). Then applying chain rule, we have(28)
h
′
t
~
=
-
c
w
f
′
T
,
h
′
′
t
~
=
-
c
w
2
f
′
′
T
.Using these relations and renaming t
~
=
t, (27) can be rewritten as(29)
m
w
2
h
′
′
t
+
ρ
c
w
a
h
′
t
+
D
1
h
t
+
D
2
h
2
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
cos
⁡
n
t
+
B
n
sin
⁡
n
t
.Note that the last summation term in (29) is now 2
π periodic.


## 3.2. Phase Analysis of the Initial Value Differential Equation

Consider the homogeneous case for (29). Dividing throughout by m
w
2 yields(30)
h
′
′
t
+
Q
1
h
′
t
+
Q
2
h
t
+
Q
3
h
2
t
=
0
,where(31)
Q
1
=
ρ
c
a
m
w
,
Q
2
=
D
1
m
w
2
,
Q
3
=
D
2
m
w
2
.Defining a new vector,(32)
H
t
=
h
~
1
t
h
~
2
t
,where h
~
1
(
t
)
=
h
(
t
) and h
~
2
(
t
)
=
h
′
(
t
); we can now rewrite (30) into two first-order equations as follows:(33)
h
~
1
′
=
h
~
2
,
h
~
2
′
=
-
Q
1
h
~
2
-
Q
2
h
~
1
-
Q
3
h
~
1
2
.It is trivial to note from the system above that there are two equilibrium points, namely, the origin (
0,0
) and the point -
Q
2
/
Q
3
,
0.Considering the first equilibrium point, that is(
0,0
), system (33) can be rewritten in a matrix form as(34)
H
′
t
=
A
H
t
+
g
t
,where(35)
H
′
t
=
h
~
1
′
h
~
2
′
,
A
=
D
h
~
0
,
0
=
0
1
-
Q
2
-
Q
1
,
h
~
=
h
~
2
-
Q
1
h
~
2
-
Q
2
h
~
1
-
Q
3
h
~
1
2
,
H
t
=
h
~
1
h
~
2
,
g
t
=
0
-
Q
3
h
~
1
2
t
.The linearized part of (34) is(36)
H
′
t
=
A
H
t
.The eigenvalues of the matrix A are two negative real numbers that can be calculated using(37)
λ
1,2
=
-
Q
1
±
Q
1
2
-
4
Q
2
2and since none of the eigenvalues of A has zero real part, the equilibrium point is hyperbolic. Hence, the behavior of nonlinear system (34) is topologically equivalent to the behavior of linearized system (36), and the approximate solution for this linear system will be good first approximation for nonlinear system (34) [20]. Note that we will restrict the values of ρ,  c,  a,  m,  w,  K
E,  K
C,  A
E,  A
C,  L, and r in order to have two distinct real eigenvalues. This means that(38)
Q
1
2
>
4
Q
2
.This condition as we will see later has biological significance in relation to the displacement of the outer wall. Analyzing the behavior of linear system, it can be seen that the origin (
0,0
) is a stable node (sink) for linear system (36). Thus, it is a stable node for nonlinear system (34) as well [20]. Therefore, in area around the origin the solution is supposed to approach the equilibrium point (see Figure 2).Figure 2
Phase portrait of the system.In order to study the behavior of the nonlinear solution near the second equilibrium point, that is(
-
Q
2
/
Q
3
,
0
), this equilibrium point has to be translated to the origin using the following transformation:(39)
h
1
^
=
h
~
1
+
Q
2
Q
3
,
h
2
^
=
h
~
2
.Using above transformation, the new system can be rewritten as(40)
h
1
^
′
=
h
2
^
,
h
2
^
′
=
-
Q
1
h
2
^
+
Q
2
h
1
^
-
Q
3
h
1
^
2
.This can be written in matrix form as(41)
H
t
^
′
=
A
^
H
t
^
+
g
t
^
,where(42)
H
t
^
′
=
h
1
^
′
h
2
^
′
,
A
^
=
D
h
^
0
,
0
=
0
1
Q
2
-
Q
1
,
h
^
=
h
2
^
-
Q
1
h
2
^
+
Q
2
h
1
^
-
Q
3
h
1
^
2
,
H
t
^
=
h
1
^
h
2
^
,
g
t
^
=
0
-
Q
3
h
1
^
2
t
.The linearized part of (41) is(43)
H
t
^
′
=
A
^
H
t
^
.The eigenvalues of the matrix A
^ are two distinct real numbers with opposite signs and can be calculated using(44)
λ
1,2
=
-
Q
1
±
Q
1
2
+
4
Q
2
2
.Therefore, the behavior of nonlinear system (41) is topologically equivalent to the behavior of linearized system (43). We have found that this equilibrium point -
Q
2
/
Q
3
,
0 is unstable saddle for the linear system; therefore, it is also an unstable saddle for the nonlinear system [20].


## 3.3. Existence of the Solution to Perturbed Nonlinear Nonhomogeneous System

Let us consider again the following second-order nonhomogeneous weakly nonlinear ODE:(45)
h
′
′
t
+
Q
1
h
′
t
+
Q
2
h
t
+
ϵ
Q
3
h
2
t
=
F
t
,where(46)
F
t
=
Q
4
a
P
m
+
D
3
+
a
Q
4
∑
n
=
1
N
A
n
cos
⁡
n
t
+
B
n
sin
⁡
n
tis 2
π periodic Fourier series and Q
4
=
1
/
m
w
2. This system can be written in vector form as(47)
h
~
1
′
=
h
~
2
,
h
~
2
′
=
-
Q
1
h
~
2
-
Q
2
h
~
1
+
F
t
-
ϵ
Q
3
h
~
1
2
,where h
~
1
(
t
)
=
h
(
t
) and h
~
2
(
t
)
=
h
′
(
t
).The associated reduced system whenϵ
=
0 is(48)
h
~
1
′
=
h
~
2
,
h
~
2
′
=
-
Q
1
h
~
2
-
Q
2
h
~
1
+
F
t
.A particular solution of reduced system (48) is expressed as(49)
ψ
t
=
Q
4
a
P
m
+
D
3
Q
2
0
+
∑
n
=
1
N
1
Q
2
-
n
2
2
+
n
Q
1
2
a
n
∗
n
b
n
∗
cos
⁡
n
t
+
b
n
∗
-
n
a
n
∗
sin
⁡
n
t
,where(50)
a
n
∗
=
Q
2
-
n
2
Q
4
A
n
-
n
Q
1
Q
4
B
n
,
b
n
∗
=
n
Q
1
Q
4
A
n
+
Q
2
-
n
2
Q
4
B
n
.Note that, from here, we can obtain all the general solutions for (48) which will require the application of the initial conditions h
~
1
(
0
)
=
h
~
2
(
0
)
=
0 to solve for any unknown constants in the general solution. Then, (45) has a 2
π periodic solution h
=
h
(
t
,
ϵ
) that is analytic in ϵ. Therefore, the solution can be expressed as convergent power series [21]:(51)
h
t
,
ϵ
=
h
0
t
+
ϵ
h
1
t
+
ϵ
2
h
2
t
+
⋯
=
∑
n
=
0
∞
ϵ
n
h
n
t
,where each function h
n
(
t
) is 2
π periodic in t.


## 3.4. Perturbation Technique to Find First-Order Approximation

Let us consider (24) with perturbed term:(52)
m
c
2
f
′
′
-
c
t
-
ρ
c
2
a
f
′
-
c
t
+
D
1
f
-
c
t
+
ϵ
D
2
f
2
-
c
t
=
F
twith the following initial condition: f
(
0
)
=
0, f
′
(
0
)
=
0, where(53)
F
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
cos
⁡
n
w
t
+
B
n
sin
⁡
n
w
tand ϵ here is assumed to be very small. Following the ideas in Lemma 1, the above differential equation can be rewritten as(54)
m
w
2
h
′
′
t
+
ρ
c
w
a
h
′
t
+
D
1
h
t
+
ϵ
D
2
h
2
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
cos
⁡
n
t
+
B
n
sin
⁡
n
twith h
(
0
)
=
0
,
h
′
(
0
)
=
0. Let the solution of differential equation (54) be expressed as power series as in the previous section. Replacing the function h
(
t
) by its power series (51) and comparing coefficients of each order, we obtain the following [22–24].Zeroth-Order Approximation. Consider(55)
m
w
2
h
0
′
′
t
+
ρ
c
w
a
h
0
′
t
+
D
1
h
0
t
=
a
P
m
+
D
3
+
a
∑
n
=
1
N
A
n
c
o
s
n
t
+
B
n
s
i
n
n
t
,
h
0
0
=
0
,
h
0
′
0
=
0
.First-Order Approximation. Consider(56)
m
w
2
h
1
′
′
t
+
ρ
c
w
a
h
1
′
t
+
D
1
h
1
t
+
D
2
h
0
t
2
=
0
,
h
1
0
=
0
,
h
1
′
0
=
0
.We can now solve (55) that yields(57)
h
0
t
=
B
~
1
e
λ
1
t
+
B
~
2
e
λ
2
t
+
a
0
+
∑
n
=
1
N
a
n
cos
⁡
n
t
+
b
n
sin
⁡
n
t
,where (58)
λ
1,2
=
-
ρ
c
w
a
±
ρ
c
w
a
2
-
4
m
w
2
D
1
2
m
w
2
,
a
0
=
a
P
m
+
D
3
D
1
,
a
n
=
a
D
1
-
n
2
w
2
m
A
n
-
a
n
ρ
w
a
c
B
n
D
1
-
n
2
w
2
m
2
+
n
ρ
w
a
c
2
,
b
n
=
a
n
ρ
w
a
c
A
n
+
a
D
1
-
n
2
w
2
m
B
n
D
1
-
n
2
w
2
m
2
+
n
ρ
w
a
c
2
,
B
~
1
=
λ
2
a
0
+
∑
n
=
0
N
a
n
-
∑
n
=
1
N
n
b
n
λ
1
-
λ
2
,
B
~
2
=
λ
1
a
0
+
∑
n
=
0
N
a
n
-
∑
n
=
1
N
n
b
n
λ
2
-
λ
1
.Since we have restriction condition from Section 3.2, that is, Q
1
2
>
4
Q
2, this is identical to the following condition:(59)
ρ
c
w
a
2
>
4
m
w
2
D
1
.Hence, λ
1
≠
λ
2
.Using solutionh
0
(
t
) from (57), we can rewrite (56) in the form(60)
m
w
2
h
1
′
′
t
+
ρ
c
w
a
h
1
′
t
+
D
1
h
1
t
+
2
D
2
B
~
1
∑
n
=
0
N
a
n
cos
⁡
n
t
e
λ
1
t
+
2
D
2
B
~
1
∑
n
=
1
N
b
n
sin
⁡
n
t
e
λ
1
t
+
2
D
2
B
~
2
∑
n
=
0
N
a
n
cos
⁡
n
t
e
λ
2
t
+
2
D
2
B
~
2
∑
n
=
1
N
b
n
sin
⁡
n
t
e
λ
2
t
+
D
2
B
~
1
2
e
2
λ
1
t
+
2
D
2
B
~
1
B
~
2
e
λ
1
+
λ
2
t
+
D
2
B
~
2
2
e
2
λ
2
t
+
D
2
a
0
2
+
2
D
2
a
0
∑
n
=
1
N
a
n
cos
⁡
n
t
+
b
n
sin
⁡
n
t
+
D
2
2
∑
n
=
1
N
a
n
2
+
D
2
2
∑
n
=
1
N
b
n
2
+
D
2
2
∑
n
=
0
N
a
n
2
cos
⁡
2
n
t
+
D
2
∑
n
=
1
N
a
n
b
n
sin
⁡
2
n
t
-
D
2
2
∑
n
=
0
N
b
n
2
cos
⁡
2
n
t
=
0
.Solution to (60) can be derived to be(61)
h
1
t
=
X
1
+
S
0
e
λ
1
t
+
X
2
+
W
0
e
λ
2
t
+
U
0
+
E
0
+
∑
n
=
1
N
U
n
cos
⁡
n
t
+
∑
n
=
1
N
Q
n
sin
⁡
n
t
+
Z
1
e
2
λ
1
t
+
Z
2
e
2
λ
2
t
+
Z
3
e
λ
1
+
λ
2
t
+
∑
n
=
1
N
E
n
cos
⁡
2
n
t
+
∑
n
=
1
N
G
n
sin
⁡
2
n
t
+
∑
n
=
1
N
S
n
cos
⁡
n
t
e
λ
1
t
+
∑
n
=
1
N
V
n
sin
⁡
n
t
e
λ
1
t
+
∑
n
=
1
N
W
n
cos
⁡
n
t
e
λ
2
t
+
∑
n
=
1
N
Y
n
sin
⁡
n
t
e
λ
2
t
,where (62)
U
0
=
-
D
2
a
0
+
∑
n
=
1
N
a
n
2
/
2
+
∑
n
=
1
N
b
n
2
/
2
D
1
,
S
0
=
-
2
D
2
B
~
1
a
0
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
,
E
0
=
-
D
2
a
0
2
2
D
1
,
W
0
=
-
2
D
2
B
~
2
a
0
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
,
U
n
=
-
2
D
2
a
0
D
1
-
n
2
m
w
2
a
n
+
ρ
c
a
n
w
b
n
D
1
-
n
2
m
w
2
2
+
ρ
c
a
n
w
2
,
Q
n
=
2
D
2
a
0
ρ
c
a
n
w
a
n
-
D
1
-
n
2
m
w
2
b
n
D
1
-
n
2
m
w
2
2
+
ρ
c
a
n
w
2
,
Z
1
=
-
D
2
B
~
1
2
4
m
w
2
λ
1
2
+
2
ρ
c
w
a
λ
1
+
D
1
,
Z
2
=
-
D
2
B
~
2
2
4
m
w
2
λ
2
2
+
2
ρ
c
w
a
λ
1
+
D
1
,
Z
3
=
-
2
D
2
B
~
1
B
~
2
m
w
2
λ
1
+
λ
2
2
+
ρ
c
w
a
λ
1
+
λ
2
+
D
1
,
E
n
=
D
2
D
1
-
4
m
n
2
w
2
/
2
b
n
2
-
a
n
2
+
2
ρ
c
a
n
w
a
n
b
n
D
1
-
4
m
n
2
w
2
2
+
2
ρ
c
a
n
w
2
,
G
n
=
-
D
2
2
ρ
c
a
n
w
/
2
a
n
2
-
b
n
2
+
D
1
-
4
m
n
2
w
2
a
n
b
n
D
1
-
4
m
n
2
w
2
2
+
2
ρ
c
a
n
w
2
,
S
n
=
2
D
2
B
~
1
2
n
w
2
m
λ
1
+
n
ρ
w
a
c
b
n
-
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
-
m
n
2
w
2
a
n
2
n
w
2
m
λ
1
+
ρ
c
a
n
w
2
+
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
-
m
n
2
w
2
2
,
V
n
=
-
2
D
2
B
~
1
2
n
w
2
m
λ
1
+
n
ρ
w
c
a
a
n
+
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
-
m
n
2
w
2
b
n
2
n
w
2
m
λ
1
+
ρ
c
a
n
w
2
+
m
w
2
λ
1
2
+
ρ
c
w
a
λ
1
+
D
1
-
m
n
2
w
2
2
,
W
n
=
2
D
2
B
~
2
2
n
w
2
m
λ
2
+
ρ
n
w
m
c
b
n
-
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
-
m
n
2
w
2
a
n
2
n
w
2
m
λ
2
+
ρ
c
a
n
w
2
+
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
-
m
n
2
w
2
2
,
Y
n
=
-
2
D
2
B
~
2
2
n
w
m
λ
2
+
n
ρ
w
m
c
a
n
+
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
-
m
n
2
w
2
b
n
2
n
w
2
m
λ
2
+
ρ
c
a
n
w
2
+
m
w
2
λ
2
2
+
ρ
c
w
a
λ
2
+
D
1
-
m
n
2
w
2
2
,
X
1
=
P
1
+
Φ
1
λ
1
-
λ
2
,
X
2
=
P
2
+
Φ
2
λ
2
-
λ
1
,where (63)
P
1
=
λ
2
-
2
λ
1
Z
1
-
λ
2
Z
2
-
λ
1
Z
3
+
λ
2
U
0
,
Φ
1
=
λ
2
∑
n
=
1
N
U
n
-
∑
n
=
1
N
2
n
G
n
+
λ
2
-
λ
1
∑
n
=
0
N
S
n
+
λ
2
∑
n
=
0
N
E
n
-
∑
n
=
1
N
n
V
n
+
Y
n
+
Q
n
,
P
2
=
-
λ
1
Z
1
+
λ
1
-
2
λ
2
Z
2
-
λ
2
Z
3
+
λ
1
U
0
,
Φ
2
=
λ
1
∑
n
=
1
N
U
n
-
∑
n
=
1
N
2
n
G
n
+
λ
1
-
λ
2
∑
n
=
0
N
W
n
+
λ
1
∑
n
=
0
N
E
n
-
∑
n
=
1
N
n
V
n
+
Y
n
+
Q
n
.For simplicity, we restricted the solution to be first-order approximation, so(64)
h
t
=
h
0
t
+
ϵ
h
1
tand then we can have a first-order approximation of the solution we are looking for which is given by(65)
u
0
,
t
=
f
-
c
t
=
f
T
=
h
t
=
h
0
t
+
ϵ
h
1
t
,where h
0 and  h
1 are defined in (57) and (61), respectively.


## 4. Linearization of the Growth and Remodeling Model

Since the growth and remodeling model discussed in the previous section is a nonlinear model that does not admit an exact solution, it would be be beneficial to consider the linearized version of this model and solve it analytically using Laplace transform. The linearization can be done using Taylor series expansion of the nonlinear terms in the boundary condition in (13) to yield(66)
u
t
t
=
c
2
u
x
x
,
(67)
u
x
,
0
=
u
t
x
,
0
=
0
,
(68)
m
u
t
t
0
,
t
=
a
P
B
t
-
K
E
A
E
L
u
0
,
t
-
K
C
A
C
L
r
2
u
0
,
t
-
K
C
A
C
1
-
r
2
2
r
2
+
ρ
c
2
a
u
x
0
,
t
,
(69)
u
t
L
,
t
=
-
c
u
x
L
,
t
.To obtain the exact solution to this linearized system, we will consider the Laplace transform of PDE (66) which becomes(70)
s
2
U
x
,
s
=
c
2
U
x
x
x
,
s
. The solution to (70) can be easily shown to be(71)
U
x
,
s
=
C
1
cosh
⁡
s
x
c
2
+
m
s
+
C
2
sinh
⁡
s
x
c
2
+
m
s
.Taking the Laplace transform of boundary conditions (68)-(69), we obtain(72)
m
s
2
U
0
,
s
=
a
P
B
-
K
E
A
E
L
U
0
,
s
-
K
C
A
C
L
r
2
U
0
,
s
-
K
C
A
C
1
-
r
2
2
r
2
+
ρ
c
2
a
U
x
0
,
s
,
(73)
s
U
L
,
s
=
-
c
U
x
L
,
s
.In order to find the constants C
1 and C
2, we substitute (71) into (73) and compare the coefficient which yields(74)
C
2
=
-
C
1
.Since we are interested in solution to displacement of cerebral spinal fluid (CSF) at x
=
0, we take Laplace transform of boundary condition at x
=
0 (72) and then substitute the solution (71) into it to obtain(75)
U
0
,
s
=
a
P
m
/
s
-
K
C
A
C
i
-
r
2
/
2
r
2
s
+
∑
n
=
1
N
a
A
n
s
/
s
2
+
n
2
w
2
+
a
B
n
n
w
/
s
2
+
n
2
w
2
m
s
2
+
ρ
c
a
s
+
K
E
A
E
/
L
+
K
C
A
C
/
L
r
2
.Taking the inverse of Laplace transformation will give us the solution of the displacement of the CSF at x
=
0 which is given by(76)
u
0
,
t
=
A
1
+
A
2
e
r
1
t
+
A
3
e
r
2
t
+
∑
n
=
1
N
C
n
e
r
1
t
+
D
n
e
r
2
t
+
E
n
cos
⁡
n
w
t
+
F
n
n
w
sin
⁡
n
w
t
,where (77)
r
1,2
=
-
ρ
c
a
±
ρ
2
c
2
a
2
-
4
m
K
E
A
E
/
L
+
K
C
A
C
/
L
r
2
2
m
,
A
1
=
a
P
m
-
K
C
A
C
1
-
r
2
/
2
r
2
m
r
1
r
2
,
A
2
=
a
P
m
-
K
C
A
C
1
-
r
2
/
2
r
2
m
r
1
r
1
-
r
2
,
A
3
=
a
P
m
-
K
C
A
C
1
-
r
2
/
2
r
2
m
r
2
r
2
-
r
1
,
C
n
=
a
A
n
r
1
+
a
B
n
n
w
m
r
1
-
r
2
r
1
2
+
n
2
w
2
,
D
n
=
a
A
n
r
2
+
a
B
n
n
w
m
r
2
-
r
1
r
2
2
+
n
2
w
2
,
E
n
=
-
C
n
-
D
n
,
F
n
=
r
2
C
n
+
r
1
D
n
+
r
1
+
r
2
E
n
.Equation (76) above represents a solution to the linearized model for the movement of the outer wall for all time t
≥
0. Obviously, the solution is the summation of exponential terms, constant terms, and periodic terms. The form of the solution suggests that the solution should start from the origin since the given initial condition is u
(
0,0
)
=
f
(
0
)
=
0. It is well known that the constant and periodic terms are bounded. Hence, the exponential term should be analyzed in order to understand the behavior of solution (76).The exponential terms in above solution (76) are in the forms e
λ
t, and note that, in order for r
1,2 to be real and distinct, we will restrict the values of ρ,  c,  a,  m,  K
E,  K
C,  A
E,  A
C,  L, and r. Hence, we need(78)
ρ
c
a
2
>
4
m
K
E
A
E
L
+
K
C
A
C
L
r
2
.Furthermore, from the given data for the problem and since all the constants in r
1,2 are positive, we found that both values of r
1,2 are negative. Therefore, the exponential terms will act in the first few time steps only and they will go to zero as t becomes large. In this case, the periodic terms and the constant will be dominant as t approaches infinity. This decay may be considered biologically significant in the displacement of the outer wall which relates to the rupture of the aneurysm.


## 5. Computational Experiments

In this section we perform numerical studies that will help to validate the theory presented in this work. For our computational experiments, the following realistic values are utilized. For the CSF,ρ
=
1000 kg/m3 and c
=
1500 m/s are used. For the wall, a
=
0.01 m2, K
E
=
800 N/m, K
C
=
3.52 N/m, A
E
=
20 m2, A
C
=
10 m2, r
=
2 m, m
=
0.001 kg, and L
=
1.5 m are used. Finally, P
m
=
8759.279403 mmHg and w
=
1 rad/s are used for the blood pressure model, and, for the harmonics, A
1
=
-
7.13, A
2
=
-
3.08, A
3
=
-
0.130, A
4
=
-
0.205, A
5
=
0.0662, B
1
=
4.64, B
2
=
-
1.18, B
3
=
-
0.564,  and B
4
=
-
0.346, B
5
=
-
0.120, all in mmHg [11]. Note that these values meet both restrictions (59) and (78), respectively. For example, (78) can be verified as follows:(79)
ρ
c
a
2
=
1000
k
g
/
m
3
∗
1500
m
/
s
∗
0.01
m
2
2
=
2.25
∗
1
0
8
k
g
2
/
s
2
>
4
m
K
E
A
E
L
+
K
C
A
C
L
r
2
=
4
∗
0.001
kg
∗
800
N
∗
20
m
2
1.5
m
+
3.52
N
∗
10
m
2
1.5
m
∗
2
2
m
=
42.6901
k
g
2
/
s
2
.Similarly (59) can also be verified.


### 5.1. Validation of the Perturbation Solution

First, we compare the first-order solution obtained in (65) using perturbation analysis against the linearized growth and remodeling model solution obtained analytically using Laplace transform in (76). In Figure 3, we plot these two solutions against the solution obtained from the linear model previously obtained by [11].Figure 3
Validation of the perturbation solution for the growth and remodeling model.First, it is interesting to note that the solution to the linearized problem with growth and remodeling obtained via Laplace transform is close to the solution for nonlinear growth and remodeling problem obtained via perturbation method. The figure also illustrates that solution for the nonlinear problem obtained via the perturbation approach also has the same behavior but differs from the solution from [11]. This goes to show that the growth and remodeling effects are important and have a significant influence on the model. Moreover, the growth and remodeling model has the effect of reducing the displacement of the outer wall and preventing the aneurysm rupture consequently. In order to understand this further, we investigate the influence of elastin and collagen contributions in the growth and remodeling process on the potential of rupture. This is discussed next.


### 5.2. Influence of Collagen and Elastin Parameters

The elastic and collagen parameters (K
E
,
K
C) seem to play an important role in the modeling of the arterial wall since they are responsible for the elasticity and strength of wall tissue. It is known from the literature that a deficiency in elastin and an increase in collagen are often related to the rupture of an aneurysm. Figure 4 shows the solution of the nonlinear coupled problem using the perturbation solution for increasing values of K
E starting from 800 N/M till 300 N/M. This corresponds to simulating the effect of deficiency in elastin and the figure clearly illustrates that as the K
E value decreases the displacement of the outer wall increases which can lead to rupture. Figure 5 shows that the displacement of the outer wall increases in a steady periodic motion as K
C increases which suggests biologically that an increase in collagen may lead to rupture. These findings are consistent with experimental findings which suggest the robustness of the one-dimensional model considered in this paper.Figure 4
Influence of elastin scaling coefficient.Figure 5
Influence of collagen scaling coefficient.


## 5.1. Validation of the Perturbation Solution

First, we compare the first-order solution obtained in (65) using perturbation analysis against the linearized growth and remodeling model solution obtained analytically using Laplace transform in (76). In Figure 3, we plot these two solutions against the solution obtained from the linear model previously obtained by [11].Figure 3
Validation of the perturbation solution for the growth and remodeling model.First, it is interesting to note that the solution to the linearized problem with growth and remodeling obtained via Laplace transform is close to the solution for nonlinear growth and remodeling problem obtained via perturbation method. The figure also illustrates that solution for the nonlinear problem obtained via the perturbation approach also has the same behavior but differs from the solution from [11]. This goes to show that the growth and remodeling effects are important and have a significant influence on the model. Moreover, the growth and remodeling model has the effect of reducing the displacement of the outer wall and preventing the aneurysm rupture consequently. In order to understand this further, we investigate the influence of elastin and collagen contributions in the growth and remodeling process on the potential of rupture. This is discussed next.


## 5.2. Influence of Collagen and Elastin Parameters

The elastic and collagen parameters (K
E
,
K
C) seem to play an important role in the modeling of the arterial wall since they are responsible for the elasticity and strength of wall tissue. It is known from the literature that a deficiency in elastin and an increase in collagen are often related to the rupture of an aneurysm. Figure 4 shows the solution of the nonlinear coupled problem using the perturbation solution for increasing values of K
E starting from 800 N/M till 300 N/M. This corresponds to simulating the effect of deficiency in elastin and the figure clearly illustrates that as the K
E value decreases the displacement of the outer wall increases which can lead to rupture. Figure 5 shows that the displacement of the outer wall increases in a steady periodic motion as K
C increases which suggests biologically that an increase in collagen may lead to rupture. These findings are consistent with experimental findings which suggest the robustness of the one-dimensional model considered in this paper.Figure 4
Influence of elastin scaling coefficient.Figure 5
Influence of collagen scaling coefficient.


## 6. Conclusion and Future Work

In this work a new mathematical model was considered that couples the elastodynamics of the blood coupled with the arterial wall dynamics and the mechanics of cerebral spinal fluid. This geometrically nonlinear coupled fluid-structure problem for the first time considered the effects of growth and remodeling and a phase analysis of the model was presented. Since this nonlinear system does not admit an exact solution, the model was analyzed using perturbation techniques to obtain approximate solutions. Section3 discusses the nonlinear analysis and perturbation solution for the model. The nonlinear coupled system was also linearized and a Laplace transform method was applied to derive the exact solution to the linearized problem. The models proposed in this paper compares well with similar models published for the linear case. Our results also seem to justify what has been experimentally observed by other research groups regarding the rupture of an aneurysm in relation to a deficiency in elastin and/or an increase in collagen. The novelty of the current paper is in the formulation, analysis, and simulation of the nonlinear coupled fluid-structure interaction model that incorporates growth and remodeling. The results from the computational experiments presented give encouragement to the authors to incorporate more features into the existing model such as incorporating viscoelasticity of the wall, adding nonlinear effects in the fluid model, and obtaining a fully implicit numerical approximation to the solution. These aspects will be considered in forthcoming papers.



---

*Source: 2869083-2016-06-16.xml*