Development of Striae Distensae - A Mechanochemical Model
Except where reference is made to the work of others, the work described in this thesis is
my own or was done in collaboration with my advisory committee. This thesis does not
include proprietary or classifled information.
Aneesh Shankar Hariharan
Certiflcate of Approval:
Amnon J Meir
Professor
Mathematics
Anotida Madzvamuse, Chair
Assistant Professor
Mathematics
Paul G Schmidt
Professor
Mathematics
Narendra Kumar Govil
Professor
Mathematics
Jo Heath
Professor
Mathematics
George.T.Flowers
Interim Dean
Graduate School
Development of Striae Distensae - A Mechanochemical Model
Aneesh Shankar Hariharan
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulflllment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
May 10, 2007
Development of Striae Distensae - A Mechanochemical Model
Aneesh Shankar Hariharan
Permission is granted to Auburn University to make copies of this thesis at its
discretion, upon the request of individuals or institutions and at
their expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Aneesh Shankar Hariharan was born on the 26th of October, 1983 in Chennai, India.
He completed his schooling at Kalpakkam, Tamilnadu, India in 2001 and then went on to
pursue his undergraduate education in Computer Science and Engineering at Pondicherry
Engineering College, Pondicherry, India. He was conferred the Bachelor of Technology
Degree in May 2005 and joined the graduate program at Auburn University, AL, USA, in
the Department of Mathematics and Statistics as a Graduate Teaching Assistant in August
2005.
iv
Thesis Abstract
Development of Striae Distensae - A Mechanochemical Model
Aneesh Shankar Hariharan
Master of Science, May 10, 2007
(B.Tech, Pondicherry Engineering College, Pondicherry, India)
101 Typed Pages
Directed by Anotida Madzvamuse
Striae distensae, otherwise known as stretch marks, are common skin lesions found
in a variety of clinical settings. They occur frequently during adolescence or pregnancy
when there is rapid tissue expansion, but may also occur in severe weight loss and in corti-
costeroid excess. Despite a considerable volume of investigative research, the pathogenesis
of striae distensae remains obscure. The interpretation of histologic samples, the major
investigative tool, demonstrates an association between dermal lymphocytic in ammation,
elastolysis and a scarring response. In this thesis, we investigate the pathogenesis of striae
distensae by addressing the coupling between mechanical forces and dermal pathology. We
develop a mathematical model for illustrating the efiect of flbroblasts and the extracellu-
lar matrix in the formation of striae, use linear stability analytical theory to predict the
range of parameter values where instability sets in and conflrm them using flnite difierence
schemes. We show how a mathematical approach provides a realistic framework that may
account for the formation of the stretch marks.
Keywords: Striae Distensae, Mechanochemical model, Stretch marks, Difiusion-driven in-
stability, Finite Difierence.
v
Acknowledgments
I would like to express my immense gratitude to Dr Anotida Madzvamuse for serving
as a role model and constantly encouraging and supporting me throughout this work. He
has been instrumental in fostering my interest in the fleld of Applied Mathematics. I
would also like to thank Professor Amnon Meir for being an inspiration and patiently going
through the initial versions of the manuscript and suggesting valuable changes. I am very
grateful to all my committee members for consenting to mentor this work. I sincerely
thank the Department of Mathematics and Statistics, Auburn University for accepting me
into its graduate program with a teaching assistantship. My father, Dr Y Hariharan and
my mother, Mrs Parimala Hariharan, have been my childhood inspiration for leading an
academic life, for which I am eternally grateful. I would also like to thank Dr A Bharathi,
for being a friend and mentor. Finally, I would like to thank my best friend and flance,
Sangeetha Srinivasan, who has encouraged me all through this work and whose company
has been invaluable to me.
vi
Style manual or journal used Journal of Approximation Theory (together with the style
known as \aums"). Bibliograpy follows van Leunen?s A Handbook for Scholars.
Computer software used The document preparation package TEX (speciflcally LATEX)
together with the departmental style-flle aums.sty.
vii
Table of Contents
List of Figures x
1 INTRODUCTION 1
1.1 Striae Distensae: An overview and motivation . . . . . . . . . . . . . . . . . 1
1.2 Ways of expressing the problem mathematically . . . . . . . . . . . . . . . . 2
1.3 Derivation of the mathematical models . . . . . . . . . . . . . . . . . . . . . 4
1.3.1 Fibroblast density conservation equation . . . . . . . . . . . . . . . . 5
1.3.2 Conservation equation for the ECM density . . . . . . . . . . . . . . 7
1.3.3 Fibroblast: ECM interaction equation . . . . . . . . . . . . . . . . . 8
1.4 Non-dimensionalization of the Model Equations . . . . . . . . . . . . . . . . 11
2 LINEAR STABILITY ANALYSIS 16
2.1 Mode selection and parameter values . . . . . . . . . . . . . . . . . . . . . . 22
3 FINITE DIFFERENCE SCHEMES 31
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.1.1 Derivation of flnite difierence schemes . . . . . . . . . . . . . . . . . 32
3.1.2 Equation (A): n(x;t) . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1.3 Equation (B): s(x;t) . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.1.4 Equation (C): v(x;t) . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.1.5 Updating u(x;t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.1.6 Computer implementation . . . . . . . . . . . . . . . . . . . . . . . . 44
4 RESULTS, CONCLUSIONS AND FUTURE DIRECTIONS 46
4.1 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.1.1 The 2-SBDF scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.1.2 The 2-SAFBDF scheme . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2 Conclusion and future directions . . . . . . . . . . . . . . . . . . . . . . . . 49
Bibliography 59
A Linear Stability Analysis: MATLAB 61
B The thomas algorithm 63
B.1 Description of the algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
B.1.1 The forward elimination phase . . . . . . . . . . . . . . . . . . . . . 64
B.1.2 The backward substitution phase . . . . . . . . . . . . . . . . . . . . 64
viii
C Implementation of the Finite Difference and Tri-Diagonal Solver:
Fortran 65
D Graphs of the solutions obtained: MATLAB 85
E Error Convergence: MATLAB 86
F Transient Solutions: MATLAB 88
ix
List of Figures
2.1 Dispersion relation for the positive root by flxing all the parameter values
other than ? as listed in the text with b=0.01,  = 0:01 and varying ?: (a)
? = 0:01, (b) ? = 0:1, (c) ? = 0:4, (d) ? = 0:5, (e) ? = 1, (f) ? = 16. . . . . 25
2.2 Dispersion relation for the positive root by flxing all the parameter values
other than ? as shown in the text with b = 0:12,  = 0:01 and varying ?: (a)
? = 0.01, (b) ? = 0:1, (c) ? = 0:6, (d) ? = 0:7, (e) ? = 1, (f) ? = 30, and (g)
? = 50. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.3 Dispersion relation for the positive root by flxing the parameter values other
than b with as listed in the text with ? = 0:5,  = 0:01 and taking (a)
b = 1, (b) b = 0:04, (c) and b = 0:001 to illustrate the efiect that a decrease
in b implies an increase in corticosteroid excess thereby giving rising to the
formation of striae. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.1 Solutions obtained by flxing all the parameters with ? = 0:01 and b=0.01
(a) n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for
n, s and u using the 2-SBDF scheme for equation (A). . . . . . . . . . . . . 51
4.2 Solutions obtained by flxing all the parameters with ? = 0:7 and b=0.01 (a)
n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for n, s
and u using the 2-SBDF scheme for equation (A). . . . . . . . . . . . . . . 52
4.3 Solutions obtained by flxing all the parameters with ? = 0:01 and b=0.01
(a) n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for
n, s and u using the 2-SAFBDF scheme for equation (A). . . . . . . . . . . 53
4.4 Solutions obtained by flxing all the parameters with ? = 0:7 and b=0.01 (a)
n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for n, s
and u using the 2-SAFBDF scheme for equation (A). . . . . . . . . . . . . . 54
4.5 Solutions obtained by flxing all the parameters with ? = 0:01 and b=0.12
(a) n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for
n, s and u using the 2-SAFBDF scheme for equation (A). . . . . . . . . . . 55
4.6 Solutions obtained by flxing all the parameters with ? = 0:7 and b=0.12 (a)
n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for n, s
and u using the 2-SAFBDF scheme for equation (A). . . . . . . . . . . . . . 56
x
4.7 Solutions obtained by flxing all the parameters with ? = 0:5 and b=1 (a) n,
s, u solution values, (b)-(d) error graph, (e)-(g) transient solutions for n, s
and u using the 2-SAFBDF scheme for Equation (A). . . . . . . . . . . . . 57
4.8 Solutions obtained by flxing all the parameters with ? = 0:7, b= 0.01, L=10
and m=1000 (a) n, s, u solution values, (b)-(d) error graph, (e)-(g) transient
solutions for n, s and u using the 2-SAFBDF scheme for Equation (A). . . 58
xi
Development of Striae Distensae - A Mechanochemical Model
Aneesh Shankar Hariharan
May 10, 2007
Chapter 1
INTRODUCTION
1.1 Striae Distensae: An overview and motivation
Striae distensae, commonly known as stretch marks is a skin condition that does not
cause any signiflcant medical problem but causes sizable distress to those afiected. Striae
distensae results when the skin is subjected to continuous stretching; increased stress is
placed on the connective tissue due to increased size of the various parts of the body. It oc-
curs on the abdomen and the breasts of pregnant women, on the shoulders of body builders,
in adolescence, and in individuals who are overweight. Skin elastolysis leads to excessive
mast cell degranulation with subsequent damage of collagen and elastin. Prolonged use
of oral or topical corticosteroids or Cushing syndrome (increased adrenal cortical activity)
leads to the development of striae. Histological flndings suggest that in the early stages,
in ammatory changes may predominate; edema is present in the dermis. In the later stages,
the epidermis becomes thin and  attened with loss of the rete ridges. The dermis has thin,
densely packed collagen bundles arranged in a parallel array horizontal to the epidermis at
the level of the papillary dermis. Elastic stains show breakage and retraction of the elastic
flbers in the reticular dermis. The broken elastic flbers curl at the sides of the striae to
form a distinctive pattern. Scanning electron microscopy shows extensive tangles of flne,
curled elastic flbers with a random arrangement. This arrangement is in contrast to nor-
mal skin, which has thick, elastic flbers with a regular distribution. When viewed by a
transmission electron microscope, the ultrastructure of elastic and collagen flbers in striae
is similar to that of healthy skin. The factors that lead to the development of striae are
1
poorly understood. No general consensus exists as to what causes striae. One suggestion
is that they develop as a result of stress rupture of the connective tissue framework. It has
also been suggested that they develop more easily in skin that has a high proportion of rigid
cross-linked collagen, for example, in early adult life. This is evident in the striae caused by
pregnancy, lactation, weight lifting, and other stressful activities. Increased adrenal cortical
activity has been implicated in the formation of striae, as in the case of Cushing syndrome.
Additionally, the cellular and extracellular matrix (ECM) alterations that mediate the clin-
ical phenotype of stretch marks remain poorly understood [1]. There have been several
detailed studies involving the efiects and even temporary treatments for this problem, but
the pathogenesis is still a mere speculation. It has been established, though, that the cause
for striae distensae could be predominantly mechanical and in certain cases such as the
Cushing?s syndrome, it is due to corticosteroid excess. Furthermore, it has also been estab-
lished that pattern forming processes of striae distensae have lots of similarities to those
of wound healing, for which a mechanochemical model has been successfully devised and
studied by Maini [2]. This also serves as motivation for devising a mathematical framework
for striae distensae.
1.2 Ways of expressing the problem mathematically
Histological samples suggest rupture of the connective tissue framework and the ECM
alterations to facilitate the phenotype of striae. The most common cells in the connective
tissue are the flbroblasts. A flbroblast is a type of cell that synthesizes and maintains the
ECM of many animal tissues. Fibroblasts provide a structural framework (stroma) for many
tissues. The main function of flbroblasts is to maintain the structural integrity of connective
2
tissue by continuously secreting precursors of the ECM. Fibroblasts have a branched cyto-
plasm surrounding an elliptical, speckled nucleus having 1 or 2 nucleoli. Though disjointed
and scattered when they have to cover a large space, flbroblasts when crowded often locally
align paralleled in clusters. Fibroblasts make collagen, glycosaminoglycans, reticular and
elastic flbers, and glycoproteins found in the extracellular matrix. In growing individuals,
flbroblasts are dividing and synthesizing the ECM. The composition of the extracellular
matrix determines the physical properties of connective tissues. Fibroblasts are morpho-
logically heterogeneous, with diverse appearances depending on their location and activity.
Though morphologically inconspicuous, ecotypically transplanted flbroblasts can often re-
tain positional memory of the location and tissue context where they had previously resided.
Fibroblasts can also migrate slowly over substratum as individual cells. While epithelial
cells form the lining of body structures it is flbroblasts and related connective tissues which
sculpt the bulk of an organism.
The basic mechanochemical model for striae relies on two of the following experimen-
tal results: (i) flbroblasts migrate within a tissue substratum, which is made up of flbrous
extracellular matrix [3], and (ii) these cells generate large traction forces during the for-
mative stages of stretch marks [4]. Hence the mechanism we shall develop will model the
mechanical interaction between the motile cells and the elastic substratum within which
they move. Fibroblast cells move by exerting forces on their surroundings, consisting of
ECM and the surface of other cells. They use their cellular protrusions, fllopodia, which
stretch out from the cells in all directions and pulling everything in its path. The biology of
these protrusions is available in [5]. As the flbroblasts move through the ECM, they deform
3
it by virtue of their traction forces. These deformations induce anisotropy efiects which in
turn afiect the cell motion. The efiects of these various coordinations results in spatially
organized cell aggregations in striae (refer to [13] and [14] for further details).
A continuum model for striae has three governing equations:
1. The conservation equation for the flbroblast density,
2. The conservation law governing the ECM and
3. The mechanical balance of forces between the cells and the ECM.
This thesis is organized as follows: In Section 1.3, we derive the mathematical models that
describe pattern formation in striae. These models are then non-dimensionalized; as illus-
trated in Section 1.4. In Section 2, we carry out detailed linear stability studies in order to
determine ranges for parameter values that will give rise to striae. The predictions of linear
stability theory are then verifled through flnite difierence schemes to provide approximate
solutions to the model equations. Finally, in Section 4, we present our numerical results,
conclusions and highlight future research ideas.
1.3 Derivation of the mathematical models
Let V ? R3 be a simply connected bounded volume at time t 2 (0;T), T > 0 and
S ? R2 be the surface closing the volume V. Throughout the course of this work, we
denote N(x;t) and S(x;t) as the flbroblast and ECM densities respectively at position
x 2 V at time t. U(x;t) denotes the displacement vector of the ECM, that is a material
point in the matrix initially at position x and undergoes a displacement to x+U. The
following equations are thus derived.
4
1.3.1 Fibroblast density conservation equation
Let c(x,t) be the density of any material. The general conservation equation says the
rate of change of amount of material in V equals the rate of  ow of material across S into
V plus the material created in V. Thus
d
dt
Z
V
c(x;t)dV = ?
Z
V
J ?dS +
Z
V
fdV; (1.1)
where J is the  ux of material and f represents the source of material which maybe a
function of c(x;t). Applying divergence theorem to the surface integral and assuming c(x;t)
is continuous and assuming V is continuous for all time t, the above equation becomes
Z
V
(@c@t +r?J ?f(c;x;t))dV = 0: (1.2)
Since the volume is arbitrary at all times, the integral must be 0 and so the conservation
equation for c(x;t) is
@c
@t +r?J = f(c;x;t): (1.3)
If classical difiusion is the process then the generalization, for example, is J = ?Drc [11].
We can also generalize the above to a case where we have several interacting materials,
species or chemicals. We will then have a vector ci(x;t) each difiusing with its own difiusion
coe?cient Di and interacting according to the vector source term f. For simplicity, we will
restrict our studies to a one-dimensional model. The above obtained equation represents
a reaction-difiusion system. We will use this conservation equation to model the flbroblast
density. Since N represents the flbroblast density, the general form of the conservation
5
equation is
@N
@t = ?r?J +M; (1.4)
where J is the  ux of cells, that is, the number of flbroblast cells crossing a unit area per
unit time and M is cell proliferation. For simplicity, we will assume that there is no cell
proliferation. The inclusion of cell proliferation should be trivial. The  ux can be expressed
as
J = Jd +Jc
where Jd = ?DrN represent the  ux due to Fick?s law. The difiusion term models random
motion, and the  ux term models directed orientation [11]. Biologically it is known that
flbroblasts move up ridges in the substratum and more speciflcally, experiments by Guido
and Tranquillo [11] show that, within oriented collagen gels, flbroblasts move preferentially
in the direction of collagen orientation by pulling themselves along the flbers. Thus the  ux
term is due to this tendency of flbroblasts to move in the direction of collagen flbers, so
that if there is a gradient of collagen, the flbroblasts tend to re-orient so as to move up that
gradient. With U(x;t) representing the displacement vector of the ECM, the convective
 ux contribution Jc is
Jc = N@U@t : (1.5)
Here the velocity of deformation of the matrix is @U@t and the amount of cells transported is
simply N times the velocity. The convective  ux is the most important contributor to the
flbroblast transport [11].
Cells also tend to disperse randomly in a homogenous isotropic medium. Classical dif-
fusion contributes a  ux term?DrN which models the random motion in which flbroblasts
6
respond to local variations in flbroblast density and tend to move down the density gradient
[19]. Also, in the case of striae, the motility of the flbroblasts is enhanced along the strain
lines and hence D(Ux) = D0 (1+ Ux) where  > 0 is a constant. Combining all the results
obtained above, we derive the conservation equation for flbroblast density as
Nt = (DNx)x ?(N Ut)x: (1.6)
It can be easily deduced therefore that the equation simplifles to
Nt = D0 Nxx +D0  (Nxx Ux + Nx Uxx)?Nx Ut ?N Utx: (1.7)
1.3.2 Conservation equation for the ECM density
The next governing equation that we model is the conservation equation for the ECM
density. Within the extracellular matrix, collagen takes the form of a flbrous network with
an elaborate structure including cross links, and consequently there is essentially no random
re-orientation of the collagen. Because the flbroblasts degrade and produce collagen, thus
reforming the network with collagen oriented in the direction of the flbroblasts, the equation
for ECM has an angular  ux term. For simplicity, in our system we do not allow any net
change in the amount of collagen; the flbroblasts simply remodel the existing network. As a
result there is only one term in the ECM equation which involves the density of ECM, the
interaction of ECM with flbroblast, and the gradient of interaction with flbroblasts. This
term is similar to the  ux term in the equation for the flbroblasts and it also contains an
additional term which arises because the rate of collagen remodeling depends on the density
of flbroblasts doing that re-modeling [19]. The ECM primarily moves by advection and we
7
mathematically model its conservation. The conservation equation for the matrix material
S(x;t) is
@S
@t +r?(SU) = P(N;S;U): (1.8)
where matrix  ux is taken to be mainly via convection and P(N;S;U) is the rate of secretion
of matrix by the cells [11]. As already discussed, secretion and degradation is thought to
play a role in certain situations involving flbroblast organization and in wound healing [15].
However, on the time scale of cell motions that we consider, we can neglect this efiect and
hence we can assume that P = 0. Also experiments by [16] indicate P = 0 during pattern
formation. Hence the conservation equation for striae is taken as
St = ?(SUt)x = ?Sx Ut ?SUtx: (1.9)
1.3.3 Fibroblast: ECM interaction equation
The composition of the ECM within which the flbroblasts move is a lot more complex
and moreover its constituents change as development proceeds. Its mechanical properties
have not been well studied either. In this case though, we are interested only in the mechan-
ical interaction between the flbroblasts and the matrix. Also the mechanical deformations
are very small [11], so as a reasonable approximation, we take the composite material of
flbroblasts plus matrix to be modeled as a linear, isotropic visco-elastic continuum with
stress tensor  . We assume that the traction forces generated by the cell are in mechanical
equilibrium with the elastic restoring forces developed in the matrix and any external forces
8
present. The mechanical force balance flbroblast-matrix equation is then [17],
r? +KF = 0; (1.10)
where F is the external force acting on the matrix (per unit matrix) and  is the stress
tensor. K is a measure of the anchoring force that resists the movement of the ECM
relative to the subcutis [11]. This equation, when applied to a spring loaded with a weight,
simply says that the applied force is balanced by the elastic forces from the extended spring.
Considering the stress tensor  that has contributions both from the ECM and the cells,
we have
 =  ECM + cell: (1.11)
The expression for a linear visco-elastic material as derived in [17], gives the stress strain
constitutive relation as
 ECM = [(?1 +?2)?t +E0?[1+?0]; (1.12)
where E0 = E1+? , ?0 = ?1?2?, ?1 and ?2 are the shear and bulk viscosities of the ECM,
? = Ux is the strain tensor, and E and ? are the Young?s modulus and the Poisson ratio
respectively. Denoting B = E 1??(1+?)(1?2?) and A = ?1 +?2, we have
 ECM = A?t +B?: (1.13)
Next we consider the contribution to the stress tensor from the flbroblast tractions, that is
 cell. The more flbroblasts there are, the greater the traction force. Also the traction force
is enhanced along the strain lines, hence the measure of the traction force as a function of
9
Ux is modeled as `(Ux) = ?0(1 + ?Ux) where the dimensionless parameter ? quantifles the
increase in traction along the strain lines [11]. If the fllopodia, with which the flbroblasts
attach to the ECM extend beyond their immediate neighborhood, as they probably do, we
need to include a long range difiusion efiect. Hence combining all these parameters, we get
the stress tensor to be
 cell = `SN; (1.14)
where ` is as deflned above and S is replaced by S + h0Sxx to include a long range efiect
in the cell traction term [11]. Finally we consider the body force F. The matrix material is
attached to a substratum of underlying tissue or the epidermis, like ropes. We model these
restraining forces as body forces proportional to the density of the ECM and the displace-
ment of the matrix from its unrestrained position and thus take F = ?SU. Combining all
the above the force balance equation (1.10) reduces to:
(A?t +B?+`SN) = KSU: (1.15)
This equation reduces to
K SU = (AUxt +BUx +?0 (1+?Ux)(S +h0 Sxx)N)x
= AUxxt +BUxx +?0 (SN)x +?0 ?(SN Ux)x +?0 h0 (N Sxx)x +?0 h0 ?(N Ux Sxx)x
= AUxxt +BUxx +?0 N Sx +?0 SNx +?0 ?N Sx Ux +?0 ?SNx Ux +?0 ?N SUxx
+?0 h0 Nx Sxx +?0 h0 N Sxxx +?0 h0 ?Nx Ux Sxx +?0 h0 ?N Uxx Sxx
+?0 h0 ?N Ux Sxxx: (1.16)
10
1.4 Non-dimensionalization of the Model Equations
Let us deflne the following new variables and time- and spatial-scales
n(^t;^x) = N?N; s(^t;^x) = S?S; u(^t;^x) = U?U; ^t = t?t and ^x = x?x: (1.17)
We can therefore compute the following derivatives
@t
@^t =
1
t?;
@x
@^x =
1
x?;
then it follows that
n^t = N
?
t? Nt; n^x =
N?
x? Nx; n^x^x =
N?
x?2 Nxx;
u^t = U
?
t? Ut; u^x =
U?
x? Ux; U^t^x =
U?
t?x? Utx; u^x^x =
U?
x?2 Uxx;
u^x^x^x = U
?
x?3 Uxxx; u^x^x^t =
U?
x?2t? Uxxt;
s^t = S
?
t? St; s^x =
S?
x? Sx; s^x^x =
S?
x?2 Sxx; s^x^x^x =
S?
x?3 Sxxx: (1.18)
Equation (1.7) can be written as
n^t = N
?
t? (D0 Nxx +D0  (Nxx Ux + Nx Uxx)?Nx Ut ?N Utx)
= D0t? N?Nxx + D0  t? N?Nxx Ux + D0  t? N?Nx Uxx ? 1t?N?Nx Ut ? 1t?N?N Utx
= D0x
?2
t?
N?
x?2 Nxx +
D0  x?2
t?
N?
x?2 Nxx
x?Ux
x? +
D0  x?
t?
N?
x? Nx
x?2 Uxx
x?2
? N
?
x? Nx
x?
t? Ut ?N
?N x?
t?x? Utx:
11
Deflne
t? = 1T ; x? = 1L; d = D0 TL2 ; U? = x? = 1L
where T and L are deflned time and length scales. Using the deflnitions in (1.17) - (1.18)
we obtain the non- dimensionalized equation
(A) n^t = dn^x^x +d (n^x^x u^x +n^x u^x^x)?(nu^t)^x : (1.19)
Equation (1.9) can be written as
StS?
x?
U?
t? =
SxS?
x?
U?Ut
t? ?S
?SU?Utx
t?x? : (1.20)
The above representation for Equation (1.20) is valid since all we have done is multiply
both sides by S?U?t?x? > 0. Using the deflnitions (1.17) - (1.18) we obtain the equation
(B) s^t = ?s^x u^t ?su^t^x: (1.21)
Similarly equation (1.16) is multiplied by U?x?2t? > 0 yields
K St?x?2 U?U =AU
?Uxxt
x?2t? +B
U?Uxx
x?2t? +?0
1
t?
Sx
x?
U?
x? N +?0
S
t?
Nx
x?
U?
x? +?0?
N
t?
Sx
x?
U?Nx
x?
+?0?St? Nxx? U
?Ux
x? +?0?
S
t?N
U?Uxx
x?2 +?0h0Nx
Sxx
x?2
U?
t? +?0h0
N
t?
Sxxx
x?2 U
?
+?0h0 ?t? Nxx? U
?Ux
x? Sxx +?0h0
N
t?
U?Uxx
x?2 Sxx +?0h0
N
t?
U?Ux
x?
Sxxx
x? :
12
Multiplying both sides by 1S0 where S? = 1S0 and deflning N? = 1N0 and applying deflnitions
(1.17) - (1.18) we obtain
KTL2su = AS
0
u^x^x^t + BS
0
Tu^x^x +?0N0Ts^xn+?0N0Tsn^x +?0N0?Tns^xu^x +?0N0?Tsn^xu^x
+?0N0?Tsnu^x^x +?0N0Th0x?2n^xs^x^x +?0N0Th0x?2ns^x^x^x +?0N0Th0x?2?n^xu^xs^x^x
+ ?0N0Th0x?2?nu^x^xs^x^x +?0N0Th0x?2?nu^xs^x^x^x: (1.22)
Dividing both sides by KTL2 and deflning
a = AK S
0 T L2
; b = BK S
0 L2
; ? = ?0 N0K L2; h = h0L2
we get the following equation
(C) au^x^x^t +bu^x^x +? (ns^x +sn^x +hn^x s^x^x +hns^x^x^x)
+? ?(ns^x u^x +sn^x u^x +snu^x^x +hn^x u^x s^x^x +hnu^x^x s^x^x +hnu^x s^x^x^x) = su: (1.23)
We can drop the^s above all the new variables and also set ? = t without any loss of
generality to obtain the following non-dimensionalized equations
nt = dnxx +d (nxx ux +nx uxx)?(nut)x ; (1.24)
st = ?sx ut ?sutx (1.25)
auxxt +buxx +? (nsx +snx +hnx sxx +hnsxxx)
+? ?(nsx ux +snx ux +snuxx +hnx ux sxx +hnuxx sxx +hnux sxxx) = su: (1.26)
13
These new equations would be our governing equations for the rest of the work. To close the
system, we prescribe initial conditions as small random perturbations around the uniform
steady state, if it exists. Boundary conditions are taken as zero- ux. Zero- ux boundary
conditions imply self-organizing processes. As it can be seen, there are a number of pa-
rameter values throughout these system of equations, each having a physical meaning even
though they are dimensionless. We set the following parameter values based on experimental
results obtained from our collaborators.
8
>>>
>>>>
>>>
>><
>>>
>>>>
>>>>
>:
A = 108Nscm?2; D0 = 10?9cm2 s?1; N0 = 104cm?3
S0 = 10?1gcm?3; ? = 1; h0 = 10?1cm2; K = 103N cm?4 S?10 ;
B = 1; 12; 100 N cm?2;
?0 = 0:001; 0:0075; 0:01; 0:045; 0:1; 0:2 N cm?2 N?10 S?10 ;
 = 0:01; 0:05:
(1.27)
Also L = 1 cm and T = 106:5 s. The values corresponding to the non-dimensional parame-
ters can be shown to be given by:
8
>>>>
>>>
>><
>>>
>>>>
>>:
d = 10?2:5; a = 10?0:5; ? = 1; h = 0:1
b = 0:01; 0:12; 1;
? = 0:01; 0:075; 0:1; 0:45; 1; 2;
 = 0:01; 0:05:
(1.28)
The parameters with variable values suggest difierent systems independent of each other
and are in relation to common situations encountered by biologists and dermatologists. For
example, increasing values for the parameter? correspond to an increase in the traction force
14
exerted by the flbroblasts on the ECM and decreasing b values correspond to a decrease
in Young?s modulus, which may occur in cases of corticosteroid excess. The rest of the
parameters are kept constant since their increase or decrease does not have much physical
signiflcance and the studies are assumed to in with relation to varying ?, b and  , since
they would have direct biological consequences which would help us to better understanding
the causes of striae distensae. The process of non-dimensionalization has many advantages,
the evident one being that the number of relevant parameters has signiflcantly reduced
dimensionless groupings, from 10 to 7. This would ease the task of determining the dynamics
of the system [11]. The reader is referred to [12] for a detailed record on the advantages of
non-dimensionalization. One thing to be noted here is that all the dimensionless parameters,
all of which are positive, can be divided into those associated with the flbroblast properties
such as ? and those related to the matrix properties, such as b,  etc. In particular we vary
b and ?, since these signify important physical concepts of the Young?s modulus and the
cell traction force respectively. Therefore, we will concentrate on these parameters in order
to get an idea of the restraints that both the flbroblasts and the ECM experience. We also
vary  to see if it has any efiect on these parameters and the patterns generated.
Although we have modeled an analytically formidable system, the conceptual frame-
work is pretty clear. We have not included all the efiects that might be relevant such as
cell proliferation and the efiects of matrix secretion. However, it is presumed that if this
model works out well, more complex features that dermatologists feel are important can be
incorporated. We have also made sure that all the essential features for the generation of
stretch marks have been incorporated and hence the model generated will be very realistic.
15
Chapter 2
LINEAR STABILITY ANALYSIS
In a seminal paper, Turing [21] showed mathematically that small random spatial
 uctuations in an otherwise well-mixed system of reacting and difiusing chemicals could
be driven unstable in the presence of difiusion to give rise to a spatial pre-pattern [11].
This process is now widely known as difiusion-driven instability. Turing?s discovery was
ground-breaking; not only was it mathematically counter-intuitive but also it gave rise to
explaining mathematically the development of form or structure in an organism.
Ourbiological problemisverymuchconcerned with patterndevelopmentand structure.
Using Turing?s theory in order to model spatial aspects observed during pattern formation
of striae, the equation system (1.24) - (1.26) must admit spatially inhomogeneous solutions.
Considering their complexity, it is highly unlikely that useful analytical solutions can be
found to these systems. However, it isknownthat muchof anypattern formation potential is
predicted by linear analysis about the uniform steady state solutions close to the bifurcation
points [11]. It is also known that such linear predictions are not infallible and must be backed
up by numerical simulations if flnite amplitude structures are required [11]. Therefore it is
worthwhile, from a practical standpoint, to carry out a detailed linear stability analysis of
the governing equations, not only as a flrst analytical step to indicate spatial pattern forming
potentialities but also because striae may involve solutions that efiectively fall within the
linear regime. The latter case is shown to be efiectively those from a nonlinear theory close
to bifurcation from uniformity [11]. So, even though this process is quite messy and unduly
complicated, if we wish to demonstrate the powerful pattern forming capabilities of striae,
16
the efiort is worthwhile. To investigate solution behavior we flrst examine the linearization
of equations (1.24) - (1.26) about the normalized steady state !0 = (n0;s0;u0)T = (1;1;0)T
looking for solutions
!1 = n(t;x)?n0; !2 = s(t;x)?s0; and !3 = u(t;x)?u0; (2.1)
which are proportional to e?t+ikx where k is the wave number and ? is the linear growth
factor. Here we assume zero  ux boundary conditions for the system of partial difierential
equations which correspond to the case of a closed cylindrical or spherical surface in three
dimensions. Zero- ux boundary conditions imply no external input (a case of self-organizing
process) which is the case in striae distensae. Hence we now have
n(t;x) = n0 e?t+ikx; s(t;x) = s0 e?t+ikx; and u(t;x) = u0 e?t+ikx: (2.2)
Using the above obtained forms, we compute the following derivatives:
nt = ?n(t;x); st = ?s(t;x); ut = ?u(t;x);
nx = ikn(t;x); sx = iks(t;x); ux = iku(t;x);
nxx = ?k2 n(t;x); sxx = ?k2 s(t;x); uxx = ?k2 u(t;x);
sxxx = ?ik3 s(t;x); uxxx = ?ik3 u(t;x);
utx = ik?u(t;x); uxxt = ?k2 ?u(t;x): (2.3)
17
The above derivatives are identical to the derivatives of the new variables !1, !2 and !3.
Let us consider equation (1.24) and linearise about the steady state:
nt = dnxx +d (nx ux)x ?(nut)x;
(!1 +n0)t = d(!1 +n0)xx +d ((!1 +n0)x (!3 +u0)x)x ?((!1 +n0)(!3 +u0)t)x and
!1t = d!1xx +d (!1x !3x)x ?(!1 !3t)x ?!3tx:
Substituting the deflnitions (2.2) and derivatives (2.3) and simplifying where non-linear
terms are ignored we have
(?+dk2)n0 +i?ku0 = 0: (2.4)
Similarly consider the equation (1.25). Linearizing about the steady state we have the
following
st = ?(sut)x = ?((!2 +1)(!3t))x = ?(!2!3t)x ?!3tx = !2t:
Substituting the derivatives and simplifying yields
?s0 +i?ku0 = 0: (2.5)
For equation (1.26), linearizing about the steady state, we obtain the following equation
su = auxxt +buxx +?(nsx +snx +hnxsxx +hnsxxx)+
??(snxux +nsxux +nsuxx +hnxuxsxx +hnuxxsxx +hnuxxsxxx):
18
Hence
(!2 +s0)(!3 +u0) = a!3xxt +b!3xx +?[(!1 +1)(!2 +1)x +(!2 +1)(!1 +1)x+
h(!1 +1)x(!2 +1)xx +h(!1 +1)(!2 +1)xxx]+
??[(!2 +1)(!1 +1)x!3x +(!1 +1)(!2 +1)x!3x+
(!1 +1)(!2 +1)(!3 +1)xx +h(!1 +1)x!3x(!2 +1)xx+
h(!1 +1)!3xx(!2 +1)xx +hn!3x(!2 +1)xxx]:
Ignoring non-linear terms, we have
a!3xxt +b!3xx +?[!2x +!1x +h!2xxx]+??[!3xx] = u0:
Re-arranging terms after substituting the values of the derivatives (2.3) and canceling com-
mon terms on both sides, we get
i? kn0 +i? ?k?hk3? s0 ??a?k2 +bk2 +? ?k2 +1? u0 = 0: (2.6)
Combining equations (2.4) - (2.6) in matrix form:
2
66
66
64
?+dk2 0 i?k
0 ? i?k
i? k i? ?k?hk3? ??a?k2 +bk2 +? ?k2 +1?
3
77
77
75
2
66
66
64
n0
s0
u0
3
77
77
75 =
2
66
66
64
0
0
0
3
77
77
75: (2.7)
19
The (3? 3) matrix is known as the Jacobian matrix. The characteristic equation for the
Jacobian matrix can be shown to be given by
????+dk2??a?k2 +bk2 +? ?k2 +1?+?k??k?h?k3???+dk2?+??2k2 = 0
Clearly ? = 0 is a trivial root. The quadratic equation
a(k)?2 +b(k)?+c(k) = 0 (2.8)
gives rise to the non trivial roots, where
a(k) = ak2; (2.9)
b(k) = (ad+? h)k4 + (? (? ?2)+b)k2 +1; (2.10)
c(k) = d? hk6 +(d? (? ?1) +bd)k4 +dk2: (2.11)
The roots of the characteristic equation (2.8) are given by
2a(k)?1;2 = ?b(k)?
p
b(k)2 ?4a(k)c(k): (2.12)
Spatially heterogeneous solutions of the linear system are characterized by the dispersion
relation ?(k) which has Re(?(k)) ? 0 but which exhibits a range of unstable modes with
Re(?(k)) > 0 for some k 6= 0. If k = 0, then ? = 0, the only spatially stable homoge-
neous solution. Although these are stable solutions, for our purposes, we are interested in
stable spatially inhomogeneous solutions. At the onset of difiusion-driven instability, the
20
stable homogeneous solution becomes unstable, instead the spatially inhomogeneous solu-
tion becomes stable thereby characterizing the phenotype for pattern formation. For this
to occur we require that Re(?(k)) > 0 to be satisfled for some k 6= 0. All the solutions
with these k?s are then linearly unstable and grow exponentially with time. These unstable
heterogeneous solutions will evolve into flnite amplitude spatially structured solutions [11].
The key assumption here is that these linearly unstable eigenfunctions which are growing
exponentially with time will be eventually bounded by the non linear terms of the reaction
difiusion system of equations and ultimately steady state spatially inhomogeneous solution
will emerge [11]. We would intuitively expect that if a conflned set exists for the kinetics,
the same set would also contain the solutions when difiusion is included [20]. So, part of
the analysis of this mechanism involves proving the presence of a conflned state within the
positive quadrant. As it can be seen from the dispersion relation, the only way a solution
with Re(?(k)) > 0 can exist is if b(k) < 0 or c(k) < 0 or both, from Descartes?s rule of signs.
Since the only negative terms involve ?, it becomes necessary that ? > 0. Biologically too,
this makes sense, since the cell traction forces are the only contributors to the aggregative
process in the matrix-flbroblast equation and hence ? > 0. Because of the central role of the
cell traction, we shall use ? as the bifurcation parameter, while varying other parameters
such as b and  . Experimentally, it is also known that, in vitro, the traction generated by
the cell increases with time for a limited period and hence ? would be the ideal parameter to
study bifurcation. The expressions for ?, b(k) and c(k), determine the domains in parameter
space where spatially inhomogeneous linearly unstable solution exist. They also give the
bifurcation surfaces in parameter space, that is, the surface which separates homogeneous
from inhomogeneous solutions. In general it is algebraically complicated to determine these
21
surfaces. However, because of the dimensionality of the parameter space, it would be of
little conceptual help in understanding the basic features of forming striae. We could make
such kinds of analysis relatively simple by neglecting one or more factors afiecting flbroblast
motion and ECM deformation. However, as already mentioned previously, we have as such
included only the essential factors that would form a realistic model. Simpler models are
also capable of generating spatial models but we are interested in the one that will solve the
problem of striae to a feasible extent and hence we will not venture much into analyzing the
domain surface spaces. With the polynomial complexity of b(k) and c(k) in the dispersion
relation, we can look forward to complex linear growth behavior.
2.1 Mode selection and parameter values
In order to carry out a detailed analysis of the mode selection and determination of the
parameter values, we plot the zeros of the characteristic equation (2.8). It must be observed
that equation (2.8) has two roots given the ? sign as shown in equation (2.12). Therefore it
is important to consider each case separately. However, here is a proof that shows that the
negative root is not appropriate for our case study. From previous discussions, it is evident
that we are interested in solutions where Re(?(k)) > 0, since these are the ones that exhibit
unstable modes and hence generate pattern. From Descartes?s rule of sign, we also inferred
that b(k) and c(k) cannot be simultaneously positive, since in this case there would be no
positive real roots to the dispersion relation. Hence consider the roots of the characteristic
equation,
2a(k)?1;2 = ?b(k)?
p
b(k)2 ?4a(k)c(k): (2.13)
22
Thetermsunderthesquarerootarepositive. Ifb(k)ispositive, then?b(k)?pb(k)2 ?4a(k)c(k)
is negative and Re(?(k)) < 0 which is of no use to us since it implies that only spatially
homogeneous solutions are stable. If b(k) is negative and c(k) is also negative, ?4a(k)c(k)
is positive and hence we once again arrive at the same conclusion with Re(?(k)) < 0. If
b(k) is negative and c(k) is positive, we will use the fact that the characteristic equation
yields to the dispersion relation leading to pattern formation only for the solution with the
largest Re(?(k)) ? 0 [13]. In the case when b(k) is negative and c(k) positive, the positive
root would most deflnitely lead to the largest positive value for the Re(?(k)). Hence, the
negative root does not serve the purpose for pattern formation in our case and we shall
not consider it. The polynomial corresponding to the positive root can be plotted easily
in MATLAB and one can then look at possible wave numbers k2n = n2 ?2, since the time
independent solution of the spatial eigenvalue problem is proportional to cos?n?xL ? which
satisfles zero- ux boundary conditions at x = 0 and x = 1. The eigenvalue in this case is
k = n? and 1=k = 1=(n?) is a measure of the wave like pattern and the eigenvalue k is the
wave number. On flnite domains there is a discrete set of possible wave numbers [11]. To
study instability we look for Re(?(k2n)) > 0 for k2? < k2n < k2+ where k2? and k2+ are the roots
of the polynomial equations in some flnite interval. Note that the order of the polynomials
obtained is greater than two and there are possibilities of having more than two real roots
to the polynomials. We however rely on MATLAB to consider possible roots that can be
isolated, whatever their number may be.
23
Positive root
Let us consider flrst the equation
2a(k)?1;2 = ?b(k)+
p
b(k)2 ?4a(k)c(k) = 0: (2.14)
At the outset, we flx all the parameter values other than ? as follows
d = 10?2:5;a = 10?0:5;b = 0:01; = 0:01;? = 1 and h = 0:1:
It can be shown that the critical value of ? for which Re(?(k2)) > 0 i.e the region where
instability starts occurring, is given by ? ? 0:4291, the critical bifurcation point (Appendix
A). By increasing the value of ? beyond this value we obtain the results illustrated in Table
2.1 and Figure 2.1 respectively. Also we can see that when the value of ? ? 16, we can
isolate wave numbers of the form n? but are not able to go beyond n = 1. Increasing ?
does not seem to have much efiect as the upper bound seems to saturate very close to the
value of ?. Besides ? signifles the cell traction force and it is usually less than 1 in normal
skin, around 0:01 [11], so increasing ? would have no meaning physically. We have included
very high values for ? in order to see the increase in regions of instability. In order for the
system to exhibit spatially inhomogeneous solutions, ? must be greater than 16 but this
has no physical meaning. Non-uniform solutions cannot be observed for ? ? 16.
Clearly, as expected, the isolated values of k do not exponentially blow up and heuristi-
cally we can conclude that they are bounded by the square root term. The fact to be noted
here is that a bifurcation to instability occurs at ? ? 0:4291. Increasing ? corresponds
24
(a)
0 1 2 3 4 5 6 7?6
?5
?4
?3
?2
?1
0
Real(lambda(k))
k
(b)
0 1 2 3 4 5 6 7?7
?6
?5
?4
?3
?2
?1
0
Real(lambda(k))
k
(c)
0 1 2 3 4 5 6 7?7
?6
?5
?4
?3
?2
?1
0
Real(lambda(k))
k
(d)
0 1 2 3 4 5 6 7?7
?6
?5
?4
?3
?2
?1
0
1
Real(lambda(k))
k
(e)
0 1 2 3 4 5 6 7?7
?6?5
?4?3
?2?1
01
23
Real(lambda(k))
k
(f)
0 1 2 3 4 5 6 7?10
0
10
20
30
40
50
60
70
80
Real(lambda(k))
k
Figure 2.1: Dispersion relation for the positive root by flxing all the parameter values other
than ? as listed in the text with b=0.01,  = 0:01 and varying ?: (a) ? = 0:01, (b) ? = 0:1,
(c) ? = 0:4, (d) ? = 0:5, (e) ? = 1, (f) ? = 16.
25
Value of ? Region of instability Wave number of form n? isolated
0.01 None None
0.075 None None
0.1 None None
0.4 None None
0.429 None None
0.4291 [2.18,2.19] None
1.00 [1.07,2.94] None
15.0 [0.26,3.14] None
16.0 [0.26,3.15] ?
50.0 [0.15,3.15] ?
100.0 [0.11,3.15] ?
Table 2.1: Regions of instability for the linearized system obtained for increasing values of
? and the wave numbers isolated.
to an increase in the traction force exerted by flbroblasts on the ECM, and with all other
parameters flxed as above, we flnd that the threshold value of ? ? 0:4291.
Now take b = 0:12 and keep all the parameters flxed other than ?. Repeat the same
exercise as above to flnd the threshold value of ? for onset of instability. The results are
generated and shown in Table 2.2 and Figure 2.2 respectively. The onset of instability for
these parameter values occurs at the threshold value of ? ? 0:624. The value of ? where
we can isolate the wave number ? is around ? ? 29 and hence for values of ? ? 29 non-
uniform spatially inhomogeneous solutions cannot be observed (Table 2.2 and Figure 2.2).
Increasing b corresponds to increase in Young?s modulus of flbroblasts.
Similarly for b = 1, we apply the same procedure. In Table 2.3, we vary b,  and search
for a region of ? to isolate wave numbers. It turns out that only ? can be isolated as before
(see Figure 4.1 for the plot of the dispersion relation). The notable fact here is that a change
in the value of  does not afiect the threshold value of ? to a great extent. This indicates
26
(a)
0 1 2 3 4 5 6 7?6
?5
?4
?3
?2
?1
0
Real(lambda(k))
k
(b)
0 1 2 3 4 5 6 7?6
?5
?4
?3
?2
?1
0
Real(lambda(k))
k
(c)
0 1 2 3 4 5 6 7?6
?5
?4
?3
?2
?1
0
Real(lambda(k))
k
(d)
0 1 2 3 4 5 6 7?6
?5
?4
?3
?2
?1
0
1
Real(lambda(k))
k
(e)
0 1 2 3 4 5 6 7?7
?6
?5
?4
?3
?2
?1
0
1
2
Real(lambda(k))
k
(f)
0 1 2 3 4 5 6 7?20
0
20
40
60
80
100
120
140
160
Real(lambda(k))
k
Figure 2.2: Dispersion relation for the positive root by flxing all the parameter values other
than ? as shown in the text with b = 0:12,  = 0:01 and varying ?: (a) ? = 0.01, (b)
? = 0:1, (c) ? = 0:6, (d) ? = 0:7, (e) ? = 1, (f) ? = 30, and (g) ? = 50.
27
Value of ? Region of instability Wave number of form n? isolated
0.01 None None
0.075 None None
0.1 None None
0.4 None None
0.623 None None
0.624 [1.96,2.03] None
1.00 [1.16,2.71] None
29.0 [0.19,3.14] None
30.0 [0.19,3.15] ?
50.0 [0.15,3.15] ?
100.0 [0.11,3.15] ?
Table 2.2: Regions of instability for the linearised system obtained for increasing values of
? and the wave numbers isolated.
that the factor contributing to the random re-orientation of the flbroblasts does not have
a signiflcant efiect on the formation of striae, since  is related to this factor. The values
of  chosen,  = 0:01 and  = 0:05, are typical values as mentioned in [19]. Combining all
the results obtained above, we can confldently say that increasing ? beyond the respective
threshold values leads to an increase in the traction force exerted by the flbroblasts on
the ECM and hence would produce stable spatially inhomogeneous solutions which in turn
produce patterns as observed in striae. These flndings will be validated through numerical
solutions as illustrated in Chapter 3.
Having dealt with the mechanical aspect for modeling the pathogenesis of striae, we
now try to reason out the experimental result of corticosteroid excess leading to striae. For
this, we flx ? = 0:05, since we know the onset of instability occurs around this value and
take  to be 0:01. We the vary b = 1; 0:12; 0:04; 0:01 as shown in Table 2.4 and Figure 2.3
respectively. Decreasing b from 1, 0:12, 0:05 and to 0:01 leads to a bifurcation to instability
28
Value of b Value of  Bifurcation Value of ? Wave number isolation value of ?
0.01 0.01 0.43 30
0.12 0.01 0.63 30
1 0.01 1.87 145
0.01 0.05 0.43 16
0.12 0.05 0.63 30
1 0.05 1.87 147
Table 2.3: Regions of ? when b and  are varied.
Value of b Region of instability Wave number of form n? isolated
1 None None
0.12 None None
0.04 [1.72,2.58] None
0.001 [1.68,2.74] None
Table 2.4: Regions of instability for the linearised system obtained for decreasing values of
b and the wave numbers isolated.
at b = 0:04 and larger instability at b = 0:01. The spatially homogeneous solution is the
only stable solution for b = 1 and 0:12. Hence decreasing b leads to instability with the
threshold value of b = 0:05. This implies a decrease in Young?s modulus of the flbroblast
cells, which in turn occurs during corticosteroid excess. Now, all that remains is to validate
predictions obtained from the linear stability theory through appropriate numerical schemes
in order to complete the framework for the pathogenesis of striae.
29
(a)
0 1 2 3 4 5 6 7?6
?5
?4
?3
?2
?1
0
Real(lambda(k))
k
(b)
0 1 2 3 4 5 6 7?7
?6
?5
?4
?3
?2
?1
0
1
Real(lambda(k))
k
(c)
0 1 2 3 4 5 6 7?7
?6
?5
?4
?3
?2
?1
0
1
Real(lambda(k))
k
Figure 2.3: Dispersion relation for the positive root by flxing the parameter values other
than b with as listed in the text with ? = 0:5,  = 0:01 and taking (a) b = 1, (b) b = 0:04, (c)
and b = 0:001 to illustrate the efiect that a decrease in b implies an increase in corticosteroid
excess thereby giving rising to the formation of striae.
30
Chapter 3
FINITE DIFFERENCE SCHEMES
3.1 Introduction
Having shown the existence of conflned sets within the positive quadrant and hav-
ing isolated regions of instability, and also the corresponding set of threshold values for
? thereon, we have now got an indication of what set of parameter values would lead to
difiusion-driven instability. Hence the linearly unstable eigenfunctions which are growing
exponentially with time will be eventually bounded by the non linear terms of the reaction
difiusion system of equations and we can ultimately expect to obtain stable spatially inho-
mogeneous steady state set of solutions. From Smoller?s result, since there exists a conflned
set in the absence of difiusion, it would also exist in the presence of difiusion. The aim of
this chapter is to flnd out these solutions numerically for difierent set of parameter values.
Due to the complexity of the model equations derived in (1.24)-(1.26), it is only pru-
dent to look for numerical approximate solutions. We will use flnite difierence schemes as
these are simpler to implement and are known to handle boundary conditions a lot easier
than flnite elements and spectral methods [10]. In order to e?ciently implement numerical
methods to the model equations (1.24)-(1.26), it is convenient to flrst assign v to be the
velocity
@u
@t = v (3.1)
31
and therefore the equations can be re-written conveniently as
nt = dnxx +d (nxxux +nxuxx)?(nv)x; (3.2)
st = ?sx v ?svx; (3.3)
avxx +buxx +? (nsx +snx +hnx sxx +hnsxxx)
+? ?(nsx ux +snx ux +snuxx +hnx ux sxx +huxx sxx +hux sxxx) = su: (3.4)
3.1.1 Derivation of flnite difierence schemes
To implement the flnite difierence schemes, approximate flnite difierence stencils to the
derivatives have to be derived. These are derived from Taylor series expansions as illustrated
below [10]. In all our derivations, it is enough to derive only the following expansions
!(x+?x;t) = !(x;t)+ ?x1! !x(x;t)+ ?x
2
2! !xx(x;t)+
?x3
3! !xxx(x;t)+O(?x
4) (3.5)
!(x??x;t) = !(x;t)? ?x1! !x(x;t)+ ?x
2
2! !xx(x;t)?
?x3
3! !xxx(x;t)+O(?x
4) (3.6)
!(x+2?x;t) = !(x;t)+2?x1! !x(x;t)+ (2?x)
2
2! !xx(x;t)+
(2?x)3
3! !xxx(x;t)+O(?x
4) (3.7)
!(x?2?x;t) = !(x;t)?2?x1! !x(x;t)+ (2?x)
2
2! !xx(x;t)?
(2?x)3
3! !xxx(x;t)+O(?x
4) (3.8)
that will be used to derive all the approximate flnite difierence stencils. We flrst derive the
forward, backward and centred flnite difierences
!x = !(x+?x;t)?!(x;t)?x +O(?x); (3.9)
!x = !(x;t)?!(x??x;t)?x +O(?x); (3.10)
32
!x = !(x+?x;t)?!(x??x;t)2?x +O(?x2): (3.11)
The forward and backward flnite difierences are only flrst order accurate while the centred
flnite difierence is second order accurate and hence yields a better accuracy. Similarly, we
derive the rest of the necessary flnite difierence approximations focussing on minimizing the
error growth and increasing the order of approximation. It can be shown that
[!(x+?x;t)+!(x;t)][z(x+?x;t)+z(x;t)] = 4!(x;t)z(x;t)+2?x(!zx +z!x)
+?x2(!zxx +!xzx +z!xx)+O(?x3); (3.12)
and
[!(x??x;t)+!(x;t)][z(x??x;t)+z(x;t)] = 4!(x;t)z(x;t)?2?x(!zx +z!x)
+?x2(!zxx +!xzx +z!xx)?O(?x3): (3.13)
Subtracting (3.13) from (3.12) we obtain the flnite difierence approximation to the derivative
@(!z)
@x = !zx +z!x =
[!(x+?x;t)+!(x;t)][z(x+?x;t)+z(x;t)]
4?x
?[!(x??x;t)+!(x;t)][z(x??x;t)+z(x;t)]4?x +O(?x2): (3.14)
The second derivative of ! can be easily obtained as
!xx = !(x+?x;t)?2!(x;t)+!(x??x;t)(?x)2 +O(?x2): (3.15)
33
It can be shown that a second order accurate flnite difierence approximation to the third
derivative can be derived by observing that
!(x+2?x;t)?2!(x+?x;t)+2!(x??x;t)?!(x?2?x;t)
= ( 83! ? 23! ? 23! + 83!)?x3!xxx +O(?x5); (3.16)
which when re-arranged conveniently yields
!xxx = !(x+2?x;t)?2!(x+?x;t)+2!(x??x;t)?!(x?2?x;t)2(?x)3 +O(?x2): (3.17)
Similar derivations can be obtained for the time-derivatives and in particular, it can be easily
shown that a second order flnite difierence scheme can be derived for the time-derivative of
the form
!t = 3!(x;t+?t)?4!(x;t)+!(x;t??t)2?t +O(?t2): (3.18)
Having developed these preliminary results, we are now in a position to apply these
flnite difierence approximations to our model equations. Also note that we will take only
those flnite difierence approximations with at least second order accuracy. It must be
noted that taking flnite difierence approximations of higher order accuracies leads to long
computational times due to the restrictions of the time-step ?t and this has not been
possible to implement on the machine we were working with. Equally important, we will
not use flrst order flnite difierence approximations simply because they fail to converge for
our model problems.
34
3.1.2 Equation (A): n(x;t)
Consider Equation 3.2
nt = dnxx +d (nxx ux +nx uxx)?(nv)x:
We consider two flnite difierence schemes for equation (A), both second order accurate in
time and space. The flrst scheme is the well-known IMEX scheme (implicit-explicit [9, 18]):
the second order flnite backward difierentiation formula (2-SBDF) and the second scheme
is a modifled version of the flrst scheme. We will denote this scheme: the second order
almost fully implicit flnite difierentiation formula (2-SAFBDF). Implicit-explicit (IMEX)
schemes use an implicit scheme to approximate the difiusion term and an explicit scheme
to approximate the reaction terms. Such schemes have been shown to be more stable and
more accurate than most schemes for reaction-difiusion systems [9, 18]. First order accurate
flnite difierence schemes blow-up for the choice of any time-steps and space-steps around
10?2 and 10?4 and even lower (results not shown) and hence we will not discuss them
in this thesis. These two feasible schemes will be devised, both of which do not allow
numerical error magniflcation. We discuss these two methods in detail for our particular
model equation below and compare the results obtained from the two schemes later on. For
the derivations below, it is important to introduce appropriate notations. We assume that
the one-dimensional mesh connectivity covering the interval ? = [0;1] is a uniform mesh of
n elements such that
0 = x0 < x1 < ??? < xn?1 < xn = 1
35
where for each i, i = 1;???n, xi = xi?1+?x, where ?x = 1n. An unstructured mesh can also
be used. For the time t 2 [0;T], T > 0, let ?t be a typical time-step such that tq = tq?1+?t,
with q = 1;2;??? and t0 = 0. Therefore we will denote by !qj = !(xj;tq) where the point
(xj;tq), (j = 0;1;??? ;n and q = 0;1;???) coincides with the flnite difierence grid points.
Also, the subscripts j and superscripts q indicate the space-step and time-step respectively.
For example, !(x+?x;t??t) will be denoted as !q?1j+1 = !(xj+1;tq?1).
The second order flnite backward difierentiation formula: 2-SBDF
This scheme has been referenced from [9, 18]. The 2-SBDF scheme can be shown to
be given by (see [9] for speciflc details)
3nq+1j ?4nqj +nq?1j
2?t = d
nq+1j+1 ?2nq+1j +nq+1j?1
?x2 +2F(n
q
j;u
q
j;v
q
j)?F(n
q?1
j ;u
q?1
j ;v
q?1
j )
where
F(nqj;uqj;vqj) = d 2?x3(nqj+1 ?2nqj +nqj?1)(uqj+1 ?uqj?1)
+ d 2?x3(nqj+1 ?nqj?1)(uqj+1 ?2uqj +uqj?1)
? 14?x[(nqj+1 +nqj)(vqj+1 +vqj)?(nqj +nqj?1)(vqj +vqj?1): (3.19)
Notice that this scheme involves calculating the values of the variable n(x;t) at time tq+1
given the values of the solutions at the previous two time-levels tq and tq?1, requiring us
to keep in memory, solutions at two difierent time-levels. Deflne ?1 = 2d ?t?x2. Re-arranging
36
terms, we get
??1nq+1j?1 +(3+2?1)nq+1j ??1nq+1j+1 =4nqj ?nq?1j +4?tF(nqj;uqj;vqj)
?2?tF(nq?1j ;uq?1j ;vq?1j ): (3.20)
Hence (3.20) can be written as
Ajnq+1j?1 +Bjnq+1j +Cjnq+1j+1 = Dj (3.21)
where j = 0;1;??? ;n and q = 0;1;??? and
Aj = ??1; Bj = (3+2?1); Cj = ??1;
and the right-hand side is deflned by
Dj = 4nqj ?nq?1j +2?t
?
2F(nqj;uqj;vqj)?F(nq?1j ;uq?1j ;vq?1j )
?
We impose zero- ux boundary conditions on the system. The coe?cients Aj, Bj, Cj how-
ever, are not functions of the solution values. Equation (3.21) gives rise to a tridiagonal
system with the coe?cients satisfying the diagonally dominant criteria
jBjj > jAjj+jCjj; j = 0;1;??? ;n:
Hence the system is in perfect structure to be solved e?ciently using the Thomas algorithm
[10]. Therefore when j = 0 for the flrst row of the tridiagonal matrix, C0 is replaced by
37
A0+C0 and similarly when j = m for the last row of the tridiagonal matrix Am is replaced
by Am+Cm for zero  ux boundary conditions. We next derive the second scheme based on
modifying the 2-SBDF scheme such that Aj, Bj, Cj are functions of the previous solution
values.
The second order almost fully implicit flnite backward difierentiation formula:
2-SAFBDF
The second scheme is an almost fully implicit scheme. The motive behind this is to
construct a scheme that exploits the implicitness from every single term in equation (B) but
at the same time, requiring at least second order accuracy both in time and space. This idea
was originally proposed and implemented in papers by Madzvamuse [9, 8]. The second order
scheme is used to approximate the time derivative (3.18) as shown in the previous section.
Similarly the second order spatial derivative is approximated implicitly using a second order
flnite difierence scheme. The remaining derivative terms are then approximated by flnite
difierence schemes such that only at least second or higher order accurate schemes are used.
The 2-SAFBDF is therefore derived as
3nq+1j ?4nqj +nq?1j
2?t = d(
nq+1j+1 ?2nq+1j +nq+1j?1
?x2 +d 
?uq
j+1 ?u
q
j?1
2?x
!?nq+1
j+1 ?2n
q+1
j +n
q+1
j?1
?x2
!
+d 
?uq
j+1 ?2u
q
j +u
q
j?1
?x2
!?nq+1
j+1 ?n
q+1
j?1
2?x
!
? 14?x
h
(nq+1j+1 +nq+1j )(vqj+1 +vqj)?(nq+1j +nq+1j?1)(vqj +vqj?1)
i
(3.22)
Observe that fully implicit schemes are used to approximate spatial derivatives as illus-
trated on the right-hand side. It is important, as observed from numerous experiments,
38
to flrst multiply by (?x)2 both sides. This is to balance the terms ranging from O(?x) to
O(?x)3. Experience has shown us that any other way results in a highly unbalanced system.
Therefore multiplying both sides by (?x)2 and re-arranging terms we obtain a tri-diagonal
system of the form
Ajnq+1j?1 +Bjnq+1j +Cjnq+1j+1 = Dj (3.23)
where
Aj = ?
h
d+d u
q
j+1 ?u
q
j?1
2?x ?d 
uqj+1 ?2uqj +uqj?1
2?x +
?x
4 (v
q
j +v
q
j?1)
i
;
Bj = 3?x
2
2?t +2d+d 
uqj+1 ?uqj?1
?x +
?x
4 (v
q
j+1 ?v
q
j?1);
Cj = ?
h
d+d u
q
j+1 ?u
q
j?1
2?x +d 
uqj+1 ?2uqj +uqj?1
2?x ?
?x
4 (v
q
j+1 +v
q
j)
i
;
and the right-hand side is given by
Dj = (4nqj ?nq?1j )?x
2
2?t: (3.24)
Clearly, now Aj, Bj and Cj are functions of the solution values at the previous time tq. This
is the key advantage of this scheme. The Thomas algorithm is then used, just as before,
taking into account the boundary conditions by replacing C0 by A0+C0 when j = 0and Am
by Am +Cm when j = m to give numerical approximate solutions. However, this equation
does not exist independently and is in uenced at every time step by the other two equations
and we impose the same kind of boundary conditions for other equations also. Details and
computer implementation of the Thomas algorithm are given in Appendix B.
39
3.1.3 Equation (B): s(x;t)
Let us consider equation (B) that models the ECM given by (3.3)
st = ?sxut ?sutx = ?@(sv)@x :
We use the results derived in Section 3.1.1 to approximate the derivatives as
sq+1j ?sqj
?t = ?v
q
j(
sq+1j+1 ?sq+1j
?x )?s
q+1
j (
vqj+1 ?vqj?1
2?x ) (3.25)
if vj is positive and
sq+1j ?sqj
?t = ?v
q
j(
sq+1j ?sq+1j?1
?x )?s
q+1
j (
vqj+1 ?vqj?1
2?x ) (3.26)
if vj is negative. The above schemes depend on whether vj is positive or negative and
depending on the case, an equivalent tri-diagonal system is formed. Since the v0s are
updated at every stage due to dependencies on the other equations, it is important that if
vj is positive, the solver automatically chooses the flnite difierence scheme. Re-arranging
the flnite difierence schemes (3.25) and (3.26) results in tri-diagonal systems who coe?cients
are given below depending on the sign of vj.
? If vj > 0 then
Aj = 0; Bj = 1+?2(vj+1 ?vj?1)?2?2vj; Cj = 2?2vj;
where ?2 = ?t2?x.
40
? If vj < 0 then
Aj = ?2?2vj; Bj = 1+?2(vj+1 ?vj ?1)+2?2vj; Cj = 0;
for j = 0;1;??? ;n. Note that the right-hand side is Dj = 0 for both schemes.
Hence we can use the same tri-diagonal solver as explained in the Appendix B with C0 =
?2?2v(j) when v(j) < 0 and j = 0 and Am = 2?2v(j) when v(j) > 0 and j = m to impose
the zero- ux boundary conditions.
3.1.4 Equation (C): v(x;t)
Consider the flnite difierence scheme for equation (C) given by (3.4).
avxx +buxxx +? (nsx +snx +hnx sxx +hnsxxx)
+? ?(nsx ux +snx ux +snuxx +hsnx ux sxx +hnuxx sxx +hnux sxxx) = su: (3.27)
Since we are only interested in the solution values of v(x;t) from this equation, a fully
implicit scheme is used to approximate the spatial derivative. There is no time-derivative
explicit in this equation. Therefore we use a centred flnite difierence scheme for vxx at
time tq+1 and compute the rest of the derivatives at the previous time tq. It is crucially
important also to balance the ?x terms. This can be achieved by simply multiplying by
(?x)2 throughout the scheme to obtain a tri-diagonal system whose coe?cients are
Aj = ?a; Bj = 2a; Cj = ?a;
41
and whose right-hand side is given by
Dj = b(uj+1 ?2uj +uj?1)
+ ??x4
h
(sj+1 +sj)(nj+1 +nj)?(sj +sj?1)(nj +nj?1)
i
+ ?h2?x(nj+1 ?nj?1)(sj+1 ?2sj +sj?1)
+ ?h2?xnj(sj+2 ?2sj+1 +2sj?1 ?sj?2)
+ ??8 (uj+1 ?uj?1)
h
(sj+1 +sj)(nj+1 +nj)?(sj +sj?1)(nj +nj?1)
i
+??nisj(uj+1 ?2uj +uj?1)
+ ??h4(?x)2(uj+1 ?uj?1)(nj+1 ?nj?1)(sj+1 ?2sj +sj?1)
+ ?? h(?x)2nj(uj+1 ?2uj +uj?1)(sj+1 ?2sj +sj?1)
+ ??h4(?x)2nj(uj+1 ?uj?1)(sj+2 ?2sj+1 +2sj?1 ?sj?2)
?siuj?x2:
The boundary conditions are implemented by replacing Cj by Aj +Cj when j = 0 and Aj
by Aj + Cj when j = m. Since jBjj = jAjj+jCjj, a small perturbation is added to Bj to
make the matrix diagonally dominant so that the same Thomas algorithm solver can be
used to solve each of the tri-diagonal systems. The algorithm is based on Gaussian forward
and backward elimination [10]. It is very e?cient because to solve per mesh point, it needs
3(add)+3(multiply)+2(divide) operations.
42
3.1.5 Updating u(x;t)
Once the main equations have been dealt with, we update the value of uj since the third
equation has given solution values for vj which have been deflned to simplify computations.
A flrst order flnite difierence scheme in time applied to the equation
v = @u@t
gives the updated scheme
uq+1j = uqj +?tvqj; j = 0;1;??? ;n: (3.28)
Initial and boundary conditions
We take initial conditions as small random perturbations around the uniform steady
state. At time t = 0, we take
nj = 1+ 0:50d0??100 sj = 1; vj = 0; uj = 0; with ? = random:
We prescribe zero- ux boundary conditions, for example, for u(x;t) these are of the form
u(?1;t) = u(1;t);
u(m+1;t) = u(m?1;t)
for all t > 0. Similar boundary conditions are applied to n, s and v.
43
3.1.6 Computer implementation
We describe the implementation of the algorithm to solve the model equations here.
1. We are interested in global model parameter values (in double precision) given by
8
>>>>
>>>
>><
>>>
>>>>
>>:
d = 10?2:5; a = 10?0:5; ? = 1; h = 0:1;
b = 0:01; 0:12; 1;
? = 0:01; 0:075; 0:1; 0:45; 1; 2;
 = 0:01; 0:05:
2. Global solution variables are n, s, v and u. Also dummy global variables are required
to store solution values at times tq and tq?1. We deflne these as nn, sn, vn and un.
3. We also need to deflne global numerical parameter values: dt, dx and length.
4. Results are stored in the four flles: nflnsol.dat (n(x;t)), sflnsol.dat (s(x;t)), uflnsol.dat
(u(x;t)) and x.dat (xj mesh points). We also store the solution values at each time
step in the flles: wnsol.m, wssol.m and wusol.m which correspond to solution values
n, s and u. These flles are huge arrays of size (m ? n) where m = T?t, and n = L?x
where T is the flnal time and L = 1 is the length of the computational interval. These
results are used in computing the error decay and the transient solutions.
5. The two sets of data are written to their respective opened flles, they are formatted to
the needed level of accuracy and they are closed once all the iterations are completed.
44
6. An error control is made by observing the convergence of the errors by calculating
jjejj =
rP?
uiq ?uq?1i
?2
where uqi and uq?1i are two successive numerical solutions
at times tq?1 and tq.
7. Each set of data values are plotted using MATLAB and graphs of the solution states
at the flnal time and the transient solution states are obtained. The numerical error
graph is also plotted from the flles with the iterative solutions (wnsol.m wssol.m and
wusol.m)
8. The time-step and the space-steps are then varied and some parameter values are
changed (based on the results obtained from linear stability theory) and difierent sets
of results are then presented.
45
Chapter 4
RESULTS, CONCLUSIONS AND FUTURE DIRECTIONS
4.1 Numerical results
4.1.1 The 2-SBDF scheme
We flrst compile the results obtained by solving the set of equations using the 2-SBDF
scheme for equation (A). Set ?t = 10?4 and ?x = 0:5?10?2 and flx all the parameter values
other than ? as follows
d = 10?2:5; a = 10?0:5; b = 0:01;  = 0:01; ? = 1; h = 0:1:
Let us take ? = 0:01 and 0:7. The corresponding results of the numerical simulations are
illustrated in Figures 4.1 - 4.2 respectively. For the above set of parameters, the bifurcation
value for instability for ? was approximately 0:43, as predicted by linear stability analysis.
As we can see, when ? = 0:01, we get solutions very close to the homogeneous steady state
solution with deviation of the order of 10?5. For ? = 0:7, the scheme converges to a uniform
steady state. This is consistent with the results obtained from linear stability analysis, since
no wave number is isolated for this value of ? and hence the uniform steady state is actually
the solution of the system at the flnal time. This scheme is stable and the numerical errors
and the transient solutions converge. It is also seen that this scheme blows up for values
of ? > 0:8 and hence we are not able to illustrate the existence of spatially inhomogeneous
solutions for large values of ?. We perform the same experiments with b = 0:12 and b = 1
and get similar results with the uniform steady state obtained at the flnal time for difierent
46
values of ? above and below the bifurcation values, as predicted by linear stability theory.
We also illustrate the results obtained using the 2-SAFBDF scheme next.
4.1.2 The 2-SAFBDF scheme
We test the 2-SAFBDF scheme applied to equation (A) in a similar way to the 2-SBDF
scheme. We choose two sets of model parameters values
d = 10?2:5; a = 10?0:5; b = 0:01;  = 0:01; ? = 1; h = 0:1:
and
d = 10?2:5; a = 10?0:5; b = 0:12;  = 0:01; ? = 1; h = 0:1:
for the same ?x and ?t as chosen for the 2-SBDF scheme and vary ? = 0:01 and 0:7. The
results are shown in Figures 4.3 - 4.6. As we can see, these numerical results are very
similar to the ones presented using the 2-SBDF scheme. For the flrst set of parameter
values chosen, the threshold value of ? was predicted as 0:43 and in the second case the
threshold value was predicted as 0:63. In both the cases we are able to see the homogeneous
steady state solution for ? = 0:01. The scheme also converges to the uniform steady state
when ? = 0:7 for both values of b and is consistent with the results obtained using linear
stability theory, since no wave number is isolated. It was also observed that a decrease in
?x and ?t yielded a similar pattern of solution behavior for the both of the above schemes,
only that it was more resolved. Also the schemes seemed to converge only if the ratio ?t(?x)2
was less than 1.3 (approximately) else they blew up. If the length L is increased from 1 to
47
10 and m is changed from 100 to 1000 to increase the number of space steps for the flrst
set of parameter values, we are able to see the onset of spatially inhomogeneous solutions
(flgure 4.8) especially for n, but they are not prominent enough to be classifled as flnite
amplitude structures. Nevertheless, we are able to see a tendency towards formation of
spatial patterns. Hence we conclude that increase in ? leads to an increase to instability
beyond the bifurcation values in each case which implies that the cell traction force is also
increased. This shows that the pathogenesis of striae is indeed mechanical!
We now turn our attention to the other experimental hypothesis for the pathogenesis of
striae, corticosteroid excess. We have linked corticosteroid excess to be characterized by a
decrease in Young?s modulus of the flbroblast which in turn is characterized by a decrease in
parameter b. We shall see how decreasing b responds to the numerical approximations. We
take b = 1 with ? flxed as 0:5 and the results are illustrated in flgure 4.7. Linear stability
theory predicted that the threshold value was 0:04 and hence values greater than this
should exhibit homogeneous state solutions which is seen clearly. It has also been observed
that when b = 0:01, the scheme converges to the uniform steady state at the flnal time
solution. Hence we have shown computationally that corticosteroid excess, characterized
by decrease in Young?s modulus in turn characterized by b, when decreased below the
instability bifurcation value of 0:04 yields patterns as seen in striae. We do not see spatially
inhomogeneous solutions that isolate a wave number since the code blows up for extremely
small values of b, where the wave number ? is isolated.
48
4.2 Conclusion and future directions
The problem of striae distensae was successfully studied and it was derived from a
realistic model that an increase in traction force beyond a certain threshold value leads
to the formation of stretch marks, accounting for the mechanical aspect of pathogenesis
of striae and decreasing the Young?s Modulus beyond a certain threshold value may also
leads to stretch mark pattern exhibition and this accounts for the experimental evidence
suggesting that corticosteroid excess, in certain cases, may cause striae . Detailed linear
stability analysis was applied to the proposed model and results were predicted as obtained
in Tables 1 - 4 in Chapter 2 in a certain range of parameter values. These results were
conflrmed using novel flnite difierence schemes. Linear stability analysis indicates that
pattern formation is possible in a certain parameter regime and this was validated by
numerical computations. We require that Re(?(k)) > 0 for some k > 0, but the results for
k=1 occur for very high values of ? which have no physical signiflcance since the typical
value of ? in normal humans is around 0.01. Hence even though we were not able to isolate
flnite amplitude spatially inhomogeneous structures due to blow up of the schemes, we
have been successful in showing that when k 6= 1, the uniform steady state is reached by
convergent schemes. Two difierent schemes were tested for the flrst equation, and both
yielded convergent and similar solutions. The main conclusions derived are signiflcant and
their biological interpretation may have direct consequences to a daily life problem faced
by millions of people.
These results have to be validated by experiments which is important before preventive
measures can be taken to resolve the problem of striae distensae. The parameter values
arrived at could serve as guidelines for experiments to search for threshold values around
49
this region. The same model can also be extended to incorporate other features. This
work is a stepping stone and we have been successful in predicting values of characteristic
parameters that would lead to pattern formation in striae and hopefully these will flnd their
way in fostering future research in this area.
50
(a)
0 50 100150200250?0.5
0
0.5
1
1.5
x?axis
n, s, u  solution values
 
 n
su
(b)
0 0.5 1 1.5 2x 1040
0.05
0.1
0.15
0.2
n errors
x?axis 
 n
(c)
0 0.5 1 1.5 2x 1040
0.5
1
1.5
2
2.5
3
3.5x 10?7
s errors
x?axis 
 s
(d)
0 0.5 1 1.5 2x 1040
0.5
1
1.5
2
2.5
3
3.5
4x 10?6
u errors
x?axis 
 u
(e) (f)
(g)
Figure 4.1: Solutions obtained by flxing all the parameters with ? = 0:01 and b=0.01 (a)
n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for n, s and u using
the 2-SBDF scheme for equation (A).
51
(a)
0 50 100150200250?0.5
0
0.5
1
1.5
x?axis
n, s, u  solution values
 
 n
su
(b)
0 0.5 1 1.5 2x 1040
0.005
0.01
0.015
n errors
x?axis 
 n
(c)
0 0.5 1 1.5 2x 1040
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1x 10?6
s errors
x?axis 
 s
(d)
0 0.5 1 1.5 2x 1042
4
6
8
10
12
14
16x 10?8
u errors
x?axis 
 u
(e) (f)
(g)
Figure 4.2: Solutions obtained by flxing all the parameters with ? = 0:7 and b=0.01 (a) n,
s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for n, s and u using the
2-SBDF scheme for equation (A).
52
(a)
0 50 100150200250?0.5
0
0.5
1
1.5
x?axis
n, s, u  solution values
 
 n
su
(b)
0 0.5 1 1.5 2x 1040
0.005
0.01
0.015
n errors
x?axis 
 n
(c)
0 0.5 1 1.5 2x 1040
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1x 10?6
s errors
x?axis 
 s
(d)
0 0.5 1 1.5 2x 1042
4
6
8
10
12
14
16x 10?8
u errors
x?axis 
 u
(e) (f)
(g)
Figure 4.3: Solutions obtained by flxing all the parameters with ? = 0:01 and b=0.01 (a)
n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for n, s and u using
the 2-SAFBDF scheme for equation (A).
53
(a)
0 50 100150200250?0.5
0
0.5
1
1.5
x?axis
n, s, u  solution values
 
 n
su
(b)
0 0.5 1 1.5 2x 1040
0.005
0.01
0.015
n errors
x?axis 
 n
(c)
0 0.5 1 1.5 2x 1040
1
2
3
4
5
6
7x 10?5
s errors
x?axis 
 s
(d)
0 0.5 1 1.5 2x 1040
0.2
0.4
0.6
0.8
1
1.2x 10?5
u errors
x?axis 
 u
(e) (f)
(g)
Figure 4.4: Solutions obtained by flxing all the parameters with ? = 0:7 and b=0.01 (a) n,
s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for n, s and u using the
2-SAFBDF scheme for equation (A).
54
(a)
0 50 100150200250?0.5
0
0.5
1
1.5
x?axis
n, s, u  solution values
 
 n
su
(b)
0 0.5 1 1.5 2x 1040
0.005
0.01
0.015
n errors
x?axis 
 n
(c)
0 0.5 1 1.5 2x 1040
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1x 10?6
s errors
x?axis 
 s
(d)
0 0.5 1 1.5 2x 1040
0.5
1
1.5x 10?7
u errors
x?axis 
 u
(e) (f)
(g)
Figure 4.5: Solutions obtained by flxing all the parameters with ? = 0:01 and b=0.12 (a)
n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for n, s and u using
the 2-SAFBDF scheme for equation (A).
55
(a)
0 50 100150200250?0.5
0
0.5
1
1.5
x?axis
n, s, u  solution values
 
 n
su
(b)
0 0.5 1 1.5 2x 1040
0.005
0.01
0.015
n errors
x?axis 
 n
(c)
0 0.5 1 1.5 2x 1040
1
2
3
4
5
6
7x 10?5
s errors
x?axis 
 s
(d)
0 0.5 1 1.5 2x 1040
0.2
0.4
0.6
0.8
1
1.2x 10?5
u errors
x?axis 
 u
(e) (f)
(g)
Figure 4.6: Solutions obtained by flxing all the parameters with ? = 0:7 and b=0.12 (a) n,
s, u solution values, (b)-(d) error graph, (e)-(g) transient solution for n, s and u using the
2-SAFBDF scheme for equation (A).
56
(a)
0 50 1001502002500
0.5
1
1.5
x?axis
n, s, u  solution values
 
 n
su
(b)
0 0.5 1 1.5 2x 1040
0.005
0.01
0.015
n errors
x?axis 
 n
(c)
0 0.5 1 1.5 2x 1040
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5x 10?5
s errors
x?axis 
 s
(d)
0 0.5 1 1.5 2x 1040
1
2
3
4
5
6
7
8x 10?6
u errors
x?axis 
 u
(e) (f)
(g)
Figure 4.7: Solutions obtained by flxing all the parameters with ? = 0:5 and b=1 (a) n, s,
u solution values, (b)-(d) error graph, (e)-(g) transient solutions for n, s and u using the
2-SAFBDF scheme for Equation (A).
57
(a)
0 500 1000 1500?0.5
0
0.5
1
1.5
x?axis
n, s, u  solution values
 
 n
su
(b)
0 0.5 1 1.5 2x 1040
2
4
6
8x 10?3
n errors
x?axis 
 n
(c)
0 0.5 1 1.5 2x 1040
0.2
0.4
0.6
0.8
1
1.2
1.4x 10?4
s errors
x?axis 
 s
(d)
0 0.5 1 1.5 2x 1040.6
0.8
1
1.2
1.4
1.6
1.8
2
2.2
2.4x 10?5
u errors
x?axis 
 u
(e) (f)
(g)
Figure 4.8: Solutions obtained by flxing all the parameters with ? = 0:7, b= 0.01, L=10
and m=1000 (a) n, s, u solution values, (b)-(d) error graph, (e)-(g) transient solutions for
n, s and u using the 2-SAFBDF scheme for Equation (A).
58
Bibliography
[1] Alaiti, Sameer. Striae Distensae. eMedicine, (2006).
[2] Maini, P. K. A Mechanochemical model for adult dermal wound contraction. J. of Biol.
Systems 3(4), 1021-1031, (1995).
[3] Hay, E. Cell Biology of the Extracellular Matrix, Newyork, Plenum, (1981).
[4] Harris, A. K., Stopak, D., Wild, D. Fibroblast traction as a mechanism for collegen
morphogenesis. Nature, 290, 249-251, (1981).
[5] Trinkaus, J. P. Cells Into Organs. Prentice-Hall, (1984).
[6] Viennet, C. Mechanical behavior of flbroblasts included in collagen lattices. J. Soc.
Biol., 195(4), 427-30 (2001).
[7] Viennet, C. Contractile forces generated by striae distensae flbroblasts embedded in
collagen lattices. Arch. Dermatol. Res., 297(1), 7-10, (2005).
[8] Madzvamuse, A. A modifled backward Euler scheme for advection-reaction-difiusion
systems. Mathematical Modeling of Biological Systems Volume I. A. Deutsch, L. Br-
usch, H. Byrne, G. de Vries and H.-P. Herzel (eds). Birkhuser, Boston, 191-197 (2007).
[9] Madzvamuse, A. Time-stepping schemes for moving grid flnite elements applied to
reaction-difiusion systems on flxed and growing domains. J. of Comp. Phys., 24 (1),
239-263, (2005).
[10] Morton, K. W. and Mayers, D. F. Numerical Solution of Partial Difierential Equations.
Cambridge University Press, (1994).
[11] Murray, J. D. Mathematical Biology. Springer Velag, (1989).
[12] Segel, Lee. Mathematics Applied to Continuum Mechanics. Dover Publications, (1987).
[13] Murray, J. D. A mechanical model for mesenchymal morphogenesis. J.Math.Biol., 17,
125-129, (1983)
[14] Murray, J. D., Oster, G.F. Generation of Biological Pattern and Form. IMA J. Math.
Appl. in Medicine and Biol., 1, 51-75, (1984)
[15] Murray,J. D., Maini,P.K., Tranquillo,R. T. Mechanical models for generating biological
patterns and form in development. Phys. Reports, 171, 60-84, (1988)
59
[16] Hinchclifie, J. R., Johnson, D. R. The development of Vertebrate Limb. Oxford Claren-
don Press, (1980)
[17] Landau, L., Lifshitz, E. Theory of Elasticity: 2nd end. Pergamon, (1970)
[18] Ruuth, S. Implicit-explicit methods for reaction-difiusion problems in pattern-
formation. J. Math. Biol., 34(2), 148-176, (1995).
[19] Sherrat, A. A Mathematical Model for Fibroblast and Collagen Orientation. Bullet.
for Math Biol., (1998).
[20] Smoller, J. Shock Waves and Reaction Difiusion equations. Springer, (1983).
[21] Turing, A. The chemical basis of morphogenesis. Phil. Trans. R. Soc. Lond. B, 237,
37-72, (1952).
60
Appendix A
Linear Stability Analysis: MATLAB
clear all; clear variables; clf; hold on;
axes(?FontSize?,16);
a=10^(-0.5); b=1; d=10^(-2.5); h=0.1; delta=1; tau=1.86;
gamma=0.01; r=1; j=1; l=1; g=1; f=1; count=1;
for k=0:0.01:40
alpha=a*k.^2;
beta=(a*d+tau*h)*k.^4+(tau*(delta-2)+b)*k.^2+1;
sigma=d*tau*h*k.^6+(d*tau*(delta-1)+b*d)*k.^4+d*k.^2;
H(g)=(-beta+real(sqrt((beta).^2-4*alpha*sigma)));
if(H(g)>0)
c(r)=g-1;
n(r)=c(r)/100;
if(r>1)
if((n(r-1)<pi)&&(n(r)>pi))
count=count+1;
end end
S(r)=H(g-1);
r=r+1;
S(r)=H(c(r-1)+1);
end
61
j=j+1;
g=g+1;
end
if(r==1)
fprintf(?No Wave Numbers Can be Isolated\n?)
else
S(r+1)=H(c(r-1)+2);
o=[S(2:r)];
fid=fopen(?isowave1.txt?,?wt?);
fprintf(fid,?The isolated wave numbers when TAU=%d are\n?,tau);
fprintf(fid,?%2.10f\n?,o);
fprintf(fid,?The isolated wave numbers along with the lower and
higher boundary limits are\n?);
fprintf(fid,?%2.10f\n?,S);
fprintf(fid,?The approxiamate interval in which these wave numbers
are situated is\n?);
fprintf(fid,?[%2.10f\t %2.10f\t]?,(-S(1)+S(2))/2,(S(r)-S(r+1))/2);
fclose(fid);
n(1)n(r-1)end count g=1
for k=0:0.01:40
plot(k,H(g),?k?);
g=g+1;
end
62
Appendix B
The thomas algorithm
The Thomas algorithm is a simplifled form of Gaussian elimination that can be used
to solve a tridiagonal systems of equations. A tridiagonal system is written as
aixi?1 +bixi +cixi+1 = di (B.1)
where a1 = 0 and cn = 0. Hence when the each of the coe?cients ai, bi, ci and di are
deflned for values of i ranging from 0 to m, which are the number of unknowns, Thomas
Algorithm can be used in solving these system of equations. In matrix form, these system
of equations are written as
2
66
66
66
66
66
66
4
b1 c1 0
a2 b2 c2
a3 b3 ?
? ? cn?1
0 an bn
3
77
77
77
77
77
77
5
2
66
66
66
66
66
66
4
x1
x2
?
?
xn
3
77
77
77
77
77
77
5
=
2
66
66
66
66
66
66
4
d1
d2
?
?
dn
3
77
77
77
77
77
77
5
The solution of these system of equations can be obtained in O(n) steps instead of 0(n3)
steps, as a result of normal Gaussian Elimination. Also this algorithm works only for
diagonally dominant matrices and this check is made. The algorithm is implemented in
two phases, a forward elimination phase in which the ai?s are eliminated and a backward
substitution phase where the solution is got.
63
B.1 Description of the algorithm
B.1.1 The forward elimination phase
? Preserve B and D(also called as forcing vectors) by copying them to new arrays for
i=0 to step m
? For i = 1 to m do
m = aib
i?1
bi = bi ?m?ci?1
di = dk ?mdi?1
B.1.2 The backward substitution phase
? Compute
xm = dmb
m
? For i = m-1 step down to 1 do
xk = dk ?ckxk+1b
k
64
Appendix C
Implementation of the Finite Difference and Tri-Diagonal Solver: Fortran
implicit none
! ------------------------------------
! Parameters
! ------------------------------------
integer i,m,simultime, timestep,snapshot, k,max, j, maxn
parameter (max=50000,maxn=10**6)
double precision paramA, paramB, paramD, delta, paramh, gamma, kappa
double precision phi, gary, mary, lary, pary, dary, cary
double precision tau, theta
double precision lambda1, lambda2
double precision dt, dx, dx2, dx3, length, finaltime,alpha,beta
double precision random,epsilon
double precision An(0:max), Bn(0:max), Cn(0:max)
double precision Dn(0:max), As(0:max), Bs(0:max)
double precision Cs(0:max), Ds(0:max), Av(0:max)
double precision Bv(0:max), Cv(0:max), Dv(0:max)
double precision err1(1:maxn), err2(1:maxn), err3(1:maxn)
double precision errn, errs, erru
double precision Fq(0:max),Fqmin(0:max)
65
double precision P1(0:max),P2(0:max)
double precision P3(0:max),P4(0:max)
double precision P5(0:max),P6(0:max)
double precision P7(0:max),P8(0:max)
double precision P9(0:max),P10(0:max)
double precision P11(0:max),P12(0:max)
double precision P13(0:max),P14(0:max)
double precision P21(0:max),P31(0:max)
double precision P41(0:max),P51(0:max)
double precision P61(0:max),P71(0:max)
double precision P81(0:max)
double precision TempBNew1(0:max), TempBNew2(0:max), TempBNew3(0:max)
double precision ResNew1(0:max),ResNew2(0:max),ResNew3(0:max)
double precision mult1,mult2, mult3
double precision n(-1:max), nn(-1:max), nold(-1:max)
double precision s(-2:max), sn(-2:max)
double precision v(-1:max), vn(-1:max), vold(-1:max)
double precision u(-1:max), un(-1:max), uold(-1:max)
! --------------------------------------
! data
! --------------------------------------
paramA= 10.0d0**(-0.50d0)
paramB=1.0d0
66
paramD=10.0d0**(-2.50d0)
delta=1.0d0
paramh=0.10d0
gamma=0.010d0
tau=0.50d0
dt=10.0d0**(-4d0)
length=1.0d0
m=210
snapshot=10
finaltime=2.0d0
! --------------------------------------
! The constant parameters
! --------------------------------------
dx= length/m
dx2= dx*dx
dx3= dx*dx*dx
simultime=INT(finaltime/dt)
print*, simultime
lambda1= 2.0d0*paramD*dt/dx2
lambda2=dt/(2.0d0*dx)
alpha=paramD*gamma/(2.0d0*dx3)
beta=1.0d0/(4.0d0*dx)
theta=(3.0d0*dx2/(2.0d0*dt))
67
kappa=paramD*gamma
phi=(dx/4.0d0)
gary=kappa/(2.0d0*dx)
mary=dx2/(2.0d0*dt)
lary=(tau*dx/4.0d0)
pary=(tau*paramh/(2.0d0*dx))
dary=(delta*tau/8.0d0)
cary=(delta*tau*paramh/(4.0d0*dx2))
! ---------------------------------------------
! Open Files for storing Results
! ---------------------------------------------
open(25, file=?wnsol.m?, status=?unknown?)
open(26, file=?wssol.m?, status=?unknown?)
open(27, file=?wusol.m?, status=?unknown?)
open(50, file=?nerr.dat?, status=?unknown?)
open(51, file=?serr.dat?, status=?unknown?)
open(52, file=?uerr.dat?, status=?unknown?)
write(25,8)
write(26,9)
write(27,10)
8 format(/,?N=[?)
9 format(/,?S=[?)
10 format(/,?U=[?)
68
! ----------------------------------
! define initial condition
! ----------------------------------
do 1 i=0,m
epsilon=random()
n(i)=1.0d0+(0.5d0-epsilon)/10.0d0
s(i)=1.0d0
u(i)=0.0d0
v(i)=0.0d0
nn(i)=0.0d0
sn(i)=0.0d0
vn(i)=0.0d0
un(i)=0.0d0
nold(i)=1.0d0+(0.50d0-epsilon)/10.0d0
uold(i)=0.0d0
vold(i)=v(i)
1 continue
! ------------------------------------
! Starting The Iterations
! -----------------------------------
do 50 timestep= 1,simultime
!do 50 timestep=1,1
do 25 i=0,m
69
An(i)=0.0d0
Bn(i)=0.0d0
Cn(i)=0.0d0
Dn(i)=0.0d0
As(i)=0.0d0
Bs(i)=0.0d0
Cs(i)=0.0d0
Ds(i)=0.0d0
Av(i)=0.0d0
Bv(i)=0.0d0
Cv(i)=0.0d0
Dv(i)=0.0d0
25 continue
! ---------------------------------------
! Zero Flux Boundary Conditions
! ---------------------------------------
n(-1)=n(1)
n(m+1)=n(m-1)
nold(-1)=nold(1)
nold(m+1)=nold(m-1)
70
s(-2)=s(2)
s(-1)=s(1)
s(m+1)=s(m-1)
s(m+2)=s(m-2)
u(-1)=u(1)
u(m+1)=u(m-1)
uold(-1)=uold(1)
uold(m+1)=uold(m-1)
v(-1)=v(1)
v(m+1)=v(m-1)
vold(-1)=vold(1)
vold(m+1)=vold(m-1)
! ------------------------------------------
! Coefficients A,B,C,D, Solving and Errors
! ------------------------------------------
do 2 i=0,m
P1(i)=(n(i+1)-2.0d0*n(i)+n(i-1))
P2(i)=(u(i+1)-u(i-1))
71
P3(i)=(n(i+1)-n(i-1))
P4(i)=(u(i+1)-2.0d0*u(i)+u(i-1))
P5(i)=(n(i+1)+n(i))
P6(i)=(v(i+1)+v(i))
P7(i)=(n(i)+n(i-1))
P8(i)=(v(i)+v(i-1))
P9(i)=(v(i+1)-v(i-1))
P10(i)=(s(i+1)+s(i))
P12(i)=(s(i)+s(i-1))
P13(i)=(s(i+1)-2.0d0*s(i)+s(i-1))
P14(i)=(s(i+2)-2.0d0*s(i+1)+2.0d0*s(i-1)-s(i-2))
P11(i)=(nold(i+1)-2.0d0*nold(i)+nold(i-1))
P21(i)=(uold(i+1)-uold(i-1))
P31(i)=(nold(i+1)-nold(i-1))
P41(i)=(uold(i+1)-2.0d0*uold(i)+uold(i-1))
P51(i)=(nold(i+1)+nold(i))
P61(i)=(vold(i+1)+vold(i))
P71(i)=(nold(i)+nold(i-1))
P81(i)=(vold(i)+vold(i-1))
!SAFBDF
if (i .eq. 0) then
Bn(i)=theta+2.0d0*paramD+2.0d0*gary*P2(i)+phi*P9(i)
72
Cn(i)=-paramD-gary*P2(i)+gary*P4(i)-phi*P8(i)-&
paramD-gary*P2(i)-gary*P4(i)+phi*P6(i)
elseif (i .eq. m) then
An(i)=-paramD-gary*P2(i)+gary*P4(i)-phi*P8(i)-&
paramD-gary*P2(i)-gary*P4(i)+phi*P6(i)
Bn(i)=theta+2.0d0*paramD+2.0d0*gary*P2(i)+phi*P9(i)
else
An(i)=-(paramD+gary*P2(i)-gary*P4(i)+phi*P8(i))
Bn(i)=theta+2.0d0*paramD+2.0d0*gary*P2(i)+phi*P9(i)
Cn(i)=-(paramD+gary*P2(i)+gary*P4(i)-phi*P6(i))
endif
Dn(i)=(-nold(i)+4.0d0*n(i))*mary
!$$$$$$ !SBDF
!$$$$$$ Fq(i)=alpha*P1(i)*P2(i)+alpha*P3(i)*P4(i)-&
!$$$$$$ beta*(P5(i)*P6(i)-P7(i)*P8(i))
!$$$$$$ Fqmin(i)=alpha*P11(i)*P21(i)+alpha*P31(i)*P41(i)-&
!$$$$$$ beta*(P51(i)*P61(i)-P71(i)*P81(i))
!$$$$$$ if (i .eq. 0) then
!$$$$$$ Bn(i)=3.0d0+2.0d0*lambda1
!$$$$$$ Cn(i)=-2.0d0*lambda1
73
!$$$$$$ elseif (i .eq. m) then
!$$$$$$ An(i)=-2.0d0*lambda1
!$$$$$$ Bn(i)=3.0d0+2.0d0*lambda1
!$$$$$$ else
!$$$$$$ An(i)=-lambda1
!$$$$$$ Bn(i)=3.0d0+2.0d0*lambda1
!$$$$$$ Cn(i)=-lambda1
!$$$$$$ endif
!$$$$$$ Dn(i)=4.0d0*n(i)-nold(i)+4.0d0*dt*Fq(i)-2.0d0*dt*Fqmin(i)
2 continue
!--- Preserve B and forcing vectors by copying them to new arrays
DO i = 0, m
TempBNew1(i) = Bn(i)
ResNew1(i) = Dn(i)
END DO
!--- Calculating solution to tridiagonal matrix
DO i = 1, m
mult1 = An(i)/TempBNew1(i-1)
TempBNew1(i) = TempBNew1(i) - Cn(i-1) * mult1
ResNew1(i) = ResNew1(i) - ResNew1(i-1) * mult1
END DO
74
nn(m) = ResNew1(m)/TempBNew1(m)
DO i = m-1, 0, -1
nn(i) = (ResNew1(i) - Cn(i) * nn(i+1)) / &
TempBNew1(i)
END DO
!error n
errn=0.0d0
do j=0,m
errn=errn+(nn(j)-n(j))**(2.0d0)
end do
err1(timestep)=sqrt(errn)
!update n
do i=0,m
nold(i)=n(i)
n(i)=nn(i)
nn(i)=0.0d0
end do
! Hyperbolic Equation: s
do 30 i=0,m
75
if(v(i) .le. 0.0d0)then
if (i .eq. 0) then
Bs(i)=1.0d0+lambda2*(v(i+1)-v(i-1))+2.0d0*lambda2*v(i)
Cs(i)=-2.0d0*lambda2*v(i)
else
As(i)=-2.0d0*lambda2*v(i)
Bs(i)=1.0d0+lambda2*(v(i+1)-v(i-1))+2.0d0*lambda2*v(i)
Cs(i)=0.0d0
endif
else
if (i .eq. m) then
As(i)=2.0d0*lambda2*v(i)
Bs(i)=1.0d0+lambda2*(v(i+1)-v(i-1))-2.0d0*lambda2*v(i)
else
As(i)=0.0d0
Bs(i)=1.0d0+lambda2*(v(i+1)-v(i-1))-2.0d0*lambda2*v(i)
Cs(i)=2.0d0*lambda2*v(i)
endif
endif
Ds(i)=s(i)
30 continue
!--- Preserve B and forcing vectors by copying them to new arrays
DO i = 0, m
76
TempBNew2(i) = Bs(i)
ResNew2(i) = Ds(i)
END DO
!--- Calculating solution to tridiagonal matrix
DO i = 1, m
mult2 = As(i)/TempBNew2(i-1)
TempBNew2(i) = TempBNew2(i) - Cs(i-1) * mult2
ResNew2(i) = ResNew2(i) - ResNew2(i-1) * mult2
END DO
sn(m) = ResNew2(m)/TempBNew2(m)
DO i = m-1, 0, -1
sn(i) = (ResNew2(i) - Cs(i) * sn(i+1)) / &
TempBNew2(i)
END DO
! errors for s(x,t)
errs=0.0d0
do j=0,m
errs=errs+(sn(j)-s(j))**(2.0d0)
end do
77
err2(timestep)=sqrt(errs)
!update s
do i=0,m
s(i)=sn(i)
sn(i)=0.0d0
end do
! Elliptic Equation
do 4 i=0,m
if (i .eq. 0) then
Bv(i)=2.0d0*paramA+0.00010d0
Cv(i)=-2.0d0*paramA
elseif (i .eq. m) then
Av(i)=-2.0d0*paramA
Bv(i)=2.0d0*paramA+0.00010d0
else
Av(i)=-paramA
Bv(i)=2.0d0*paramA+0.00010d0
Cv(i)=-paramA
endif
78
Dv(i)=paramB*P4(i)+lary*(P10(i)*P5(i)-P12(i)*P7(i)) &
+pary*P3(i)*P13(i)+pary*n(i)*P14(i)&
+dary*P2(i)*(P10(i)*P5(i)-P12(i)*P7(i)) &
+8.0d0*dary*n(i)*s(i)*P4(i)+cary*P2(i)*P3(i)*P13(i) &
+4.0d0*cary*n(i)*P4(i)*P13(i)+cary*n(i)*P2(i)*P14(i) &
-s(i)*u(i)*dx2
4 continue
! --- Preserve B and forcing vectors by copying them to new arrays
DO i = 0, m
TempBNew3(i) = Bv(i)
ResNew3(i) = Dv(i)
END DO
!--- Calculating solution to tridiagonal matrix
DO i = 1, m
mult3 = Av(i)/TempBNew3(i-1)
TempBNew3(i) = TempBNew3(i) - Cv(i-1) * mult3
ResNew3(i) = ResNew3(i) - ResNew3(i-1) * mult3
END DO
vn(m) = ResNew3(m)/TempBNew3(m)
79
DO i = m-1, 0, -1
vn(i) = (ResNew3(i) - Cv(i) * vn(i+1)) / &
TempBNew3(i)
END DO
!Update v
do i=0,m
vold(i)=v(i)
v(i)=vn(i)
vn(i)=0.0d0
end do
!Calculate u
do 6 i=0,m
un(i)=v(i)*dt+u(i)
6 continue
! errors for u(x,t)
erru=0.0d0
do j=0,m
erru=erru+(un(j)-u(j))**(2.0d0)
end do
err3(timestep)=sqrt(erru)
80
!Update u
do i=0,m
uold(i)=u(i)
u(i)=un(i)
un(i)=0.0d0
end do
! -------------------------------------
! Store the results in a file
! -------------------------------------
k=mod(timestep,snapshot)
if(k.eq.0)then
write(25,42) (n(i),i=0,m)
write(26,43) (s(i),i=0,m)
write(27,44) (u(i),i=0,m)
42 format(1505f15.8,?;?)
43 format(1505f15.8,?;?)
44 format(1505f15.8,?;?)
endif
write(50,77) err1(timestep)
write(51,78) err2(timestep)
81
write(52,79) err3(timestep)
77 format(f15.8)
78 format(f15.8)
79 format(f15.8)
do i=0,m
if(ABS(u(i)).gt.5.0d0) then
print*,?Blow up for u values:u(i)=?,ABS(u(i))
print*,?Number of iterations for blow up=?, timestep
stop
elseif(ABS(v(i)).gt.5.0d0) then
print*,?Blow up for v values:v(i)=?,ABS(v(i))
print*,?Number of iterations for blow up=?, timestep
stop
elseif(ABS(n(i)).gt.5.0d0) then
print*,?Blow up for n values:n(i)=?,ABS(n(i))
print*,?Number of iterations for blow up=?, timestep
stop
elseif(ABS(s(i)).gt.5.0d0) then
print*,?Blow up for s values:s(i)=?,ABS(s(i))
print*,?Number of iterations for blow up=?, timestep
82
stop
end if
end do
50 continue
close(53)
close(54)
close(55)
close(60)
close(61)
! Write the output of the solutions at the final time
open(70, file=?nfinsol.dat?, status=?unknown?)
open(71, file=?sfinsol.dat?, status=?unknown?)
open(72, file=?ufinsol.dat?, status=?unknown?)
write(70,74) (n(i),i=0,m)
write(71,75) (s(i),i=0,m)
write(72,76) (u(i),i=0,m)
74 format(f15.8)
75 format(f15.8)
76 format(f15.8)
83
close(70)
close(71)
close(72)
write(25,200)
write(26,201)
write(27,202)
200 format(/,?];?)
201 format(/,?];?)
202 format(/,?];?)
close(25)
close(26)
close(27)
stop
end
84
Appendix D
Graphs of the solutions obtained: MATLAB
clear variables; clf; close all; format long;
axes(?FontSize?,18);
hold on;
load nfinsol.dat;
load sfinsol.dat;
load ufinsol.dat;
plot(nfinsol,?r.?);
plot(sfinsol,?b:?);
plot(ufinsol,?k?);
ylabel(?n, s, u solution values?);
xlabel(?x-axis?);
85
Appendix E
Error Convergence: MATLAB
clear variables;
clf;
close all;
format long;
axes(?FontSize?,18);
load nerr.dat;
load serr.dat;
load uerr.dat;
figure(1);
plot(nerr,?r.?);
ylabel(?n errors?);
xlabel(?x-axis?);
legend(?n?);
figure(2);
plot(serr,?b?);
ylabel(?s errors?);
xlabel(?x-axis?);
legend(?s?);
figure(3);
plot(uerr,?k?);
86
ylabel(?u errors?);
xlabel(?x-axis?);
legend(?u?)
87
Appendix F
Transient Solutions: MATLAB
clear variables;
close all; clf; whitebg(?w?)
axes(?FontSize?,18);
fprintf(?READING INPUT PLOTING VALUES \n?);
fprintf(?1 is to plot n solution \n?);
fprintf(?2 is to plot s solution \n?);
fprintf(?3 is to plot u solution \n?);
nsu = input(?Type your choice \n?);
finaltime=2.0;
snapshot=10;
dt=10^(-4);
m=100;
length=1.0;
dx=length/m;
for i=0:(finaltime/(snapshot*dt))-1
time(i+1)=i*dt;
end
for j=0:m
x(j+1)=j*dx;
end
88
hold on;
fprintf(?Begin plotting \n?);
if nsu == 1
wnsol;
[r,c]=size(N);
ONE=ones(c,1);
Time=time?*ONE?;
ONE1=ones(r,1);
X=x?*ONE1?;
surfl(X?,10*Time,N);
shading interp;
colormap autumn;
elseif nsu == 2
wssol;
[r,c]=size(S);
ONE=ones(c,1);
Time=time?*ONE?;
ONE1=ones(r,1);
X=x?*ONE1?;
surfl(X?,10*Time,S)
shading interp;
colormap autumn;
elseif nsu == 3
89
wusol;
[r,c]=size(U);
ONE=ones(c,1);
Time=time?*ONE?;
ONE1=ones(r,1);
X=x?*ONE1?;
surfl(X?,10*Time,U)
shading interp;
colormap autumn;
end
ylabel(?Time t?);
xlabel(?Domain x(t)?);
90