A numerical approach to model the dynamics in a mammalian cell exposed to a carcinogenic compound

(1)

A numerical approach to model the dynamics in a mammalian cell exposed to a carcinogenic compound

Degree Project in Engineering Physics, SA104X

SARA HAMIS HAMIS@KTH.SE

School of Engineering Sciences

Royal Institute of Technology, Stockholm, Sweden Supervisor: Michael Hanke

Examiner: Mårten Olsson

April 2014

(2)

(3)

Abstract

When a mammalian cell is exposed to a carcinogen, the carcinogen gives rise to new compounds. In this paper a numerical approach is taken to find the amount of these compounds over time. The concentration of the compounds are assumed to satisfy the reaction and/or the diffusion equation, techniques used for solving these equations include the one-dimensional finite element method (FEM) and the lumped parameter approach. The mathematically derived results are compared to in-vitro measurements. The mathematical model set up in this paper will prove insufficient.

(4)

Referat

En numerisk ansats att modellera dynamiken i en däggdjusrcell som utsättsfö r en carcinogen

När en däggdjurscell utsätts för en carcinogen skapas nya äm- nen. I denna rapport tas en numerisk ansats till att modellera mängden av dessa ämnen i cellen över tid. Koncentrationen av ämnena antas satisfiera reaktions- och/eller diffusionsekva- tionen, metoder för att lösa dessa ekvationer inkluderar fini- ta element metoden (FEM) samt lumped parameter-metoden.

De matematiskt erhållna resultaten jämförs med mätningar från ett in-vitro-experiment. Den matematiska modellen framtagen i denna rapport visar sig vara otillräcklig.

(5)

Introduction

In modern medicine cancer is an important research area as it is a common and widely spread disease occurring in multiple forms. Cancer is caused when cells with damaged DNA divide and grow uncontrollably [6]. This damage can in many cases be derived from environmental factors, such as the presence of a carcinogen. Exposure to carcinogens disrupts the intracellular dynamics and may turn a healthy cell cancerous.

In a previous in-vitro experiment carried out at Karolinska Institutet in Stockholm, researchers investigated how healthy mammalian cells react when their environment is exposed to a carcinogen, namely polyclinic aromatic hydrocarbons diol epoxides (PAH DEs) [2][1]. In the in-vitro experiment, mammalian cells of the type V79 are placed in an aqueous medium into which the carcinogen is planted. The carcinogen causes several physical and chemical reactions to occur in the extra-cellular medium and in the cell. The chemical reactions give rise to new compounds including DNA-adducts, PAH tetrols and glutathione conjugates (DE GSH conjugates). DNA-adducts are cancerous while previous research [2] shows that PAH tetrols and GSH-conjugates prevent cancer as they eliminate the carcinogen. In the in-vitro experiment, the amount of three compounds (PAH DEs, PAH Tetrols and GSH conjugates) are measured respectively during a time-span of 30 minutes.

The aim of this project is to try to reenact the in-vitro experiment mathematically and thus to compute the amount of different compounds present during a time-span of 30 minutes. The concentration of the regarded compounds are assumed to satisfy the diffusion-reaction equation. When attempting to pursue a mathematical model for the in-vitro experiment, one finds it is sufficient to regard a single cell in a finite volume of extra-cellular medium due to the experimental setup where the extra-cellular volume by far exceeds the volume of all cells combined. Consequently, the mathematical model is set up in this project describes the governing equations on a closed system (one cell and a finite volume of extra-cellular medium). The equations are assembled into a system of equations and computed in the programing environment matlab. Mathematical techniques used for obtaining the system of equations include the finite element method (FEM) and the lumped parameter approach.

(8)

CHAPTER 1. INTRODUCTION

More than every third person living in Sweden today will receive a cancer diagnosis at some point during their lifetime, according to the Swedish Cancer Society, an inde- pendent non-profit organization and the largest financier of the country’s cancer research [6]. This statement alone is enough to encourage research on the subject cancer and further, computational approaches to model biological systems are increasing in popularity [8]. This increment is partly explained by advancing technology such as medical imaging and bioengineering in nano-scales, providing better insight to bio-chemical phenomena and intracellular data.

1.1 Structure of the paper

The working procedure of the project and this paper is illustrated in Figure 1.1 below.

Part I of this paper provides background-information needed to set up a mathematical model. This includes biological information regarding the cell anatomy and the in- vitro experiment. A physical background is also given, describing diffusion and chemical reactions. Finally the mathematical background introduces the FEM, and especially how to find an approximate solution to the diffusion-reaction equation using this technique. In Part II of the report the mathematical model is set up resulting in a system of equations that is computed in MATLAB. The computed results are compared to measurements from the in-vitro experiment and discussed in Part III.

Figure 1.1: Working procedure of the report.

2

(9)

Part I

Background

(10)

(11)

Chapter 2

Biological Background

This section provides the biological background information needed to mathematically model a cell in the in-vitro experiment [2].

2.1 Anatomy of the mammalian cell

The anatomy of a mammalian cell is highly complex and needs to be simplified for upcoming mathematical procedure. The DNA is confined in the cell-nucleus and as damaged DNA causes cancer, the nucleus is of special importance. Therefore the cell anatomy is simplified to consisting of two aqueous domains; the cell nucleus and the part of the cell that is not the nucleus. The latter part is called the cytoplasm and comprises various constituents such as organelles and molecules of different sizes, thus the cytoplasm is highly heterogeneous. The nucleus is protected by the nucleus membrane and the cytoplasm is protected by the cell membrane. Cells of type V79 that are used in the in-vitro experiment are assumed to have a roughly spherical shape. Now the extra cellular medium, into which the carcinogen is planted, is also a included in the mathematical model. Throughout this paper the aqueous domains are denoted i

and membranes are denoted iaccording to Table 2.1. An imaginary surface is also introduced, confining the finite volume of extra-cellular medium surrounding the cell.

The domains are illustrated in Figure 2.1.

Notation Domain in cell

Imaginary surface defining the system

1 Extracellular medium

1 Cell membrane

3 Cytoplasm

2 Nucleus membrane

5 Cell nucleus

Table 2.1: Notations for the domains constituting the system.

(12)

CHAPTER 2. BIOLOGICAL BACKGROUND

2.2 The regarded compounds

The carcinogen gives rise to new compounds, the compounds of interest are listed in Table 2.2. All compounds are denoted by a capital letter throughout the paper.

Notation Biological compound Present in domains C PAH diol epoxides 1, ₃, ₅

U PAH tetrols 1, ₃, ₅

B DE GSH conjugates 3

A DNA adducts 5

E Protein adducts 3, 5

Table 2.2: Notations for the regarded compounds.

2.3 Impact of the carcinogen in the in-vitro experiment

A cell consist of different aqueous sub-compartments separated by protecting membranes, but due to being lipophilic the carcinogen (C) diffuses through the membranes and thus reach the entire cell [2]. In all aqueous sub-compartments the carcinogen undergoes hydrolysis, ergo reacts with the water present, to produce PAH tetrols (U). The PAH tetrols are also lipophilic and hence diffuse through the membranes. What is more; in the cytoplasm the carcinogen reacts with protein present to form protein adducts (E), and with enzymes to form glutathione conjugates (B). Both latter reactions protect the cell by eliminating the carcinogen. Continuing, the carcinogen and the PAH tetrols reach the cell nucleus, containing the DNA. Here the carcinogen reacts with the DNA to forming DNA-adducts (A) which are cancerous pieces of DNA sharing covalent bonds with carcinogenic chemicals. Figure 2.1 provides a schematic view of how the compounds live in the cell.

6

(13)

2.3. IMPACT OF THE CARCINOGEN IN THE IN-VITRO EXPERIMENT

Figure 2.1: Figure a) represents the domains listed in Table 2.1. Figure b) shows the diffusion over domains (double arrows) and the chemical reactions (single arrows) occurring in the cell due to the carcinogen’s presence.

(14)

(15)

Chapter 3

Physical Background

The presence of the carcinogen causes diffusion and chemical reactions to occur, these processes are described in this section. Using bracket notation, the concentration of compound C is denoted [C], and so on for all compounds.

3.1 Diffusion

In all aqueous domains, substances U and C diffuse. Diffusion is the movement of molecules from an area of high concentration to an area of low concentrations. C and U will continue to diffuse until equilibrium is established. The general diffusion equation describes the diffusion occurring in mathematically,

ˆ[C]

ˆt = Ò(Di· Ò[C]), (3.1a)

ˆ[U]

ˆt = Ò(Di· Ò[U]), (3.1b)

on i where Di is a diffusion constant, i = 5, 3, 1 and the Ò-operator is a differen- tial operator that in Cartesian coordinates is ^û_ˆx+^û_ˆy+^û_ˆz.

3.2 Transmission

Substances C and U diffuse through the membranes 1 and 2 due to being lipophilic, diffusion through membranes is called transmission.

(16)

CHAPTER 3. PHYSICAL BACKGROUND

The transmission rate over some infinitely thin surface, S, separating 2 volumes Vaand V_b is directly proportional to the difference of concentration in Vaand Vbas illustrated in Figure 3.1. The transmission is directed from the volume with higher concentration to the volume with lower concentration.

Figure 3.1: Transmission through an infinitely thin surface S from a volume Va with high concentration to a volume Vbwith low concentration.

3.3 Chemical Reactions

The regarded chemical reactions in the in-vitro experiment are irreversible and are as follows [2],

On 1

Ó

C≠≠æ U,^r^U,1

On 3

Y_ _] __ [

C≠≠æ U,^r^U,3 C≠≠æ B,^r^B,3 C≠≠æ E,^r^E,3

On 5

Y_ _] __ [

C≠≠æ U,^rÛ,5 C≠≠æ A,^rÂ,5 C≠≠æ E,^rÊ,5

where the superscripts over the arrows denote the reaction constants.

10

(17)

3.3. CHEMICAL REACTIONS

Note that the chemical reactions cause substance C to dissolve, thus substance C has a negative reaction constants. Mathematically the reactions in the system are described by the Equations 3.2, 3.3, 3.4.

On 1,

d[U]

dt = rU,1[C], (3.2a)

d[C]

dt = ≠rU,1[C], (3.2b)

on 1,

d[U]

dt = rU,3[C], (3.3a)

d[B]

dt = rB,3[C], (3.3b)

d[E]

dt = rE,3[C], (3.3c)

d[C]

dt = ≠(rU,3+ rB,3+ rE,3)[C], (3.3d) and on 5,

d[U]

dt = rU,5[C], (3.4a)

d[A]

dt = rA,5[C], (3.4b)

d[E]

dt = rE,5[C], (3.4c)

d[C]

dt = ≠(rU,5+ rA,5+ rE,5)[C], (3.4d)

where the reaction constants have the unit [s^≠1].

(18)

(19)

Chapter 4

Mathematical Background

This section provides an introduction to the FEM and shows how the method may be used to solve partial differential equations, such as the transient diffusion-reaction equation which compounds in the in-vitro experiment obey.

4.1 The Finite Element Method

Using the FEM it is possible to simplify a complicated exact function, known or unknown, by a sum of simpler functions [3],[4],[5]. The exact function, defined on some domain , is denoted uexactand u is the the approximate function. Problems studied in this paper are one-dimensional and the independent spatial variable is denoted x, thus in the static case u may be written as a finite sum,

u(x) =^ÿ

n=1

c_nÏ_n(x), (4.1)

where Ïn(x) are some functions and cn are the respective coefficients. In this paper the functions Ïn(x) are chosen to be linear polynomials and further u is required to be continuous on , thus u is a continuous piecewise linear function as shown in Figure 4.1.

Figure 4.1: The function u consists of is a piecewise linear polynomials and acts to approximate some exact function uexact.

(20)

CHAPTER 4. MATHEMATICAL BACKGROUND

4.2 Defining the mesh

Let the domain be defined between to points and a and b on the x-axis,

= [a, b]. (4.2)

Now is divided into N non-overlapping sub-domains, or elements, so that

= ¹fi ²fi ... fi ^N, (4.3)

where ⁿ denotes element number n. In the case of linear polynomials each element is confined between two points on the x-axis called global nodes. Since there are N elements it follows that there are N + 1 global nodes,

Global nodes on : x1, x₂, ..., x_N+1,

where xn denotes global node number n. Especially x1 lies on x = a and xN+1 lies on x = b. Using nodes, the interval of an element may be written

n= [xn, xn+1]. (4.4)

Throughout this paper all nodes are set to be equidistant and separated a distance hso that

h= xn+1≠ xn. (4.5)

All concepts introduced in this subsection are demonstrated in Figure 4.2.

14

(21)

4.3. LAGRANGE POLYNOMIALS

Figure 4.2: The domain consists of N elements and N + 1 nodes, separated a distance h. On node number n, that is on xn, the function u takes the value u(xn).

4.3 Lagrange Polynomials

On every element ⁿ, u is approximated by a linear polynomial. From linear algebra one knows that a set of polynomials of maximum degree p constitute a vector space spanned by p + 1 basis vectors. Taking advantage of this idea the polynomial-functions Ïn(x) in Equation 4.5 are referred to as the basis functions, spanning a vector space V ,

V = span{Ï1, Ï2, ..., ÏN+1}. (4.6)

Now, there are different choices of such basis functions one can make and here nodal basis functions are used so that each function Ïn(x) is described by the values they take

(22)

CHAPTER 4. MATHEMATICAL BACKGROUND at the N + 1 global nodes xi. One nodal basis is the Lagrange basis spanned by the functions { ⁄1, ⁄₂, ...} with the property that ⁄i(x) equals 1 exactly on node i and ⁄i(x) equals zero on all other nodes,

⁄_i(xj) = ”ij=

I 1 for i = j , 0 otherwise.

The basis functions Ïn(x) in Equation are now set to be the Lagrange basis functions.

This means that each Ïn(x) is only non-zero on elements containing node xn and that two basis functions are non-zero on each element. The basis functions are referred to as hat-functions due to their hat-like shape demonstrated in Figure 4.3.

Figure 4.3: The hat-functions, Ïn(xi) equal 1 on node xnand 0 on all other nodes. For example Ï2equals 1 on node number 2 and 0 on all other nodes.

From Figure 4.3 one sees that the hat-functions Ïnmay be described as

16

(23)

4.4. THE LOCAL FUNCTIONS

Ï_n(x) = Y_ ] _[

x≠xⁿ≠1

h x_n≠1Æ x Æ xn,

xn+1≠x

h x_nÆ x Æ xn+1,

0 otherwise.

Using the Lagrange basis, evaluating the approximate function u on all the nodes xi

yields the following sum

u(xi) =^N+1^ÿ

n=1

cnÏn(xi), (4.7)

moreover, if x œ [xn, x_n+1] then

u(x) = cnxn+1≠ x

h + cn+1x≠ xn

h . (4.8)

4.4 The local functions

In order to simplify the upcoming calculations on each elements, the elements are now viewed locally. A new independent spatial variable denoted › is now introduced and an element is confined between 2 local nodes ›1and ›2so that

Âⁿ= [›1, ›₂] = [≠1, 1], (4.9)

for all elements ⁿ, the tilde-sign denotes that the domain is regarded locally. In the same way local functions will be denotedÏÂ_n(›). Each element is associated with 2 non- zero global functions that are translated to local functions in this section. A local basis function is simply a straight line satisfying

I ÏÂ0(›) equals 1 for › = ≠1 and 0 for › = 1, Â

Ï₁(›) equals 0 for › = ≠1 and 1 for › = 1,

as shown in Figure 4.4.

(24)

Figure 4.4: Local functions viewed on an element.

It is straightforward that the local basis functions may be written

Â Ï0(›) =1

2(1 ≠ ›), (4.10a)

Â Ï₁(›) =1

2(1 + ›). (4.10b)

4.5 Galerkin’s method

Now that the FEM mesh is defined and the basis functions are chosen it is time to show how the FEM can be used find an approximate function u to an unknown exact function u_exact. When faced with transient problems Equation may be modified to

u(x, t) =^N+1^ÿ

i=1

c_n(t)Ïn(xn), (4.11)

18

(25)

4.5. GALERKIN’S METHOD

where the time-dependence is gathered in the constants. A differential equation may be expressed using an operator L such that

L(U) = 0, (4.12)

where U is some solution to the differential equation described by the operator L. For an exact function Equation 4.12 will hold and the right-side of the equation will equal zero. However, for an approximate function u this is not the case (if the exact function is not a linear polynomial), thus the right-hand side of Equation 4.12 will not zero everywhere on the domain . The amount of which u fails to fulfill Equation 4.12 is called the residual R,

R= L(u). (4.13)

It is reasonable to state that the best approximate solution u to describe uexact is the one that minimizes the residual R. Now, the residual is a function of u and u is s a linear combination of the functions Ï1, Ï₂, ..., Ï_N+1 that span a vector-space V . Further if u satisfies some boundary-conditions (BCs) u lies in a function-space VB.C,

uœ VB.C= {B.C satisfied|u œ V } (4.14)

Using Galerkin’s method, one may minimize the residual by forcing the scalar prod- uct of the residual and some test-function v to be zero. In other words the residual and the test-function v are orthogonal. Using Galerkin’s method v is chosen to lie in the vector-space V (but not satisfying BCs) so that

v= v(x) œ V. (4.15)

The scalar product of the residual and some test-function is now set to be equal to zero,

(R, v) = 0 ’v œ V. (4.16)

(26)

This statement is easily validated by the fact that any v œ V can be expressed as a linear combination of the basis vectors spanning V ,

v=^N^ÿ⁺¹

n=1

c_nÏ_n, (4.17)

where cnœ R are some arbitrary constants, therefore on may write

(R, v) = (R,^N+1^ÿ

n=1

c_nÏ_n) =^N+1^ÿ

n=1

c_n(R, Ïn), (4.18)

which must hold for all partial sums and therefore Equation 4.18 reduces to a single term so that

(R, v) = cn(R, Ïn) = 0, (4.19)

which is true for all cnœ R (including cn”= 0) so consequently

(R, Ïn) = 0, n = 1, 2, ..., N + 1. (4.20)

The scalar product of may be written in integral form,

(R, v) =^⁄ RvdV = 0. (4.21)

4.6 The diffusion equation using the FEM

In this paper the diffusion equation and the reaction equation are of special interest as they are satisfied by the compounds in the in-vitro experiment. This section shows

20

(27)

4.6. THE DIFFUSION EQUATION USING THE FEM

the transient diffusion equation on a spherical symmetrical domain is handled using the FEM.

4.6.1 Galerkin’s method

When working with the diffusion-equation on a spherical domain L is set to be

L(U) =ˆU

ˆt ≠ Ò(ÒU) = 0), (4.22)

where all diffusion-constants are skipped for abbreviation. The residual for an approxi- mate function u becomes

R(u) =ˆu

ˆt≠ Ò(Òu)). (4.23)

Using Galerkin’s method the scalar-product of R(u) and a test-function v œ V is re- quired to equal zero so that,

(R(u), v) = 0, (4.24)

or written more elaborately in integral-form,

⁄ (ˆu(x, t)

ˆt ≠ Ò(Òu))vdV = 0. (4.25)

By moving the last term to the right-hand-side of the equation one may write

⁄ ˆu(x, t)

ˆt vdV =^⁄ Ò(Òu)vdV. (4.26)

(28)

4.6.2 Integration by Parts

The right-hand-side of Equation 4.26 is evaluated using partial integration,

⁄

Ò(Òu)vdV =

⁄

ˆ ÒuvndS ≠

⁄

ÒuÒvdV, (4.27)

where ˆ the boundary of the volume and n is the outwards facing normal of this boundary so that one may write

⁄

Ò(Òu)vdV =

⁄

ˆ

ˆu ˆndS≠

⁄

ÒuÒvdV, (4.28)

where the first term on the right-hand-side is some boundary condition. By inserting Equation 4.29 into Equation 4.26 the diffusion equation becomes an ODE,

ˆ ˆt

⁄

uvdV = ≠^⁄ ÒuÒvdV + {some B.Cs}, (4.29)

where u = u(x, t) and v = v(x) .

4.6.3 A spherical domain

If the regarded domain is of spherical kind the volume-integral over is suitably evaluated using spherical coordinates,

⁄

· dV =

⁄_b

x=a

⁄_2ﬁ

◊=0

⁄_ﬁ

„=0· x²sin(„) dx d◊ d„. (4.30)

Further if the domain is spherical symmetrical and thus only depends on the ra- dial component x the volume integral over simplifies to

⁄

· dV = 4ﬁ

⁄_b

x=a· x²dx. (4.31)

22

(29)

Now if the diffusion equation, Equation 4.29, is evaluated on a spherical symmetrical domain the nabla-operator, Ò, reduces to ^ˆu_dx so that, disregarding BCs, the equation becomes

4ﬁˆ ˆt

⁄_b

x=ax²u(x, t)v(x)dx = ≠4ﬁ^⁄ ^b

x=ax²ˆu(x, t) ˆx

ˆv(x, t)

ˆx dx, (4.32)

where a factor 4ﬁ can be divided from both sides of the equation. Remember now that the functions u and v are linear combinations of the basis functions Ï1, ..., Ï_N+1so that

u(x, t) = c(t)^N^ÿ⁺¹

n=1

Ïn(x), (4.33a)

v(x) =^N^ÿ⁺¹

m=1

Ï_n(x), (4.33b)

inserting this into the ODE-type diffusion equation, Equation 4.32 yields

N+1ÿ

n=1 N+1ÿ

m=1

d

dtcn(t)^⁄ ^b

a x²Ïn(x)Ïm(x)dx = ≠^N+1^ÿ

n=1 Nÿ+1 m=1

cn(t)^⁄^b

a x²dÏn(x) dx

dÏm(x) dx dx.

(4.34)

4.6.4 Evaluating the functions on one element using local functions The integral-equation 4.34 is now evaluated separately on each domain in accordance with the fact that

⁄

· =

⁄

1· +

⁄

2· + ... +

⁄

N

·. (4.35)

Each element will be evaluated using the local parameters from Section 4.4 so it is

(30)

now necessary to translate the global integrals in Equation into local integrals.

Locally the element ⁿ ranges from ≠1 Æ › Æ 1 and thus has an interval of length 2. Globally ⁿranges from xnÆ x Æ xn+1 and has an interval of length h, so that

dx d› =h

2, (4.36a)

d›

dx= 2

h. (4.36b)

Using Equations 4.36 an integral over an element can be transformed from global to local coordinates,

⁄

n· dx =

⁄_x_n+1

xn · dx =

⁄_›₂

›1 · dx d›d›=^⁄ ¹

≠1·v h

2d›. (4.37)

The occurring factors in the integrands in Equation 4.34 are x², the local basis functions and their derivatives. The derivatives of the basis functions are obtained by straight for- ward differentiation,

Â Ï₀(›) =1

2(1 ≠ ›), (4.38a)

Â Ï1(›) =1

2(1 + ›), (4.38b)

dÏÂ0(›)

d› = ≠1, (4.38c)

dÏÂ1

d› = 1. (4.38d)

Using the chain-rule, the derivatives with respect to x are written as derivatives with respect to ›,

dÏ(x)

dx = dÏ(›) d›

d›

dx= {using Equation 4.36 } =dÏ(›) d›

2

h. (4.39)

24

(31)

Since the x²factor occurs in the integral, x is now regarded as a function of ›,

x= x(›). (4.40)

Now, a function (here x(›)) takes on exactly the function value at a node so

x(›) =ÏÂ0(›)xn+ÏÂ1(›)xn+1, (4.41) on some domain ⁿ= [xn, xn+1]. Now local parameters are used to evaluate the integral Equation 4.34 on a single element, the equation to solve becomes

ÿ1 n=0

ÿ1

m=0˙cn(t)^⁄¹

≠1x²(›)ÏÂn(›)ÏÂm(›)h

2d›= ≠^ÿ¹

n=0

ÿ1 m=0

cn(t)^⁄¹

≠1x²(›)ÏÂ^Õ_n(›)2 hÏ^Â^Õ_m(›)2

h h 2d›, (4.42)

where ˙cn(t) denotes the time-derivative of cn(t) andÏÂ^Õ denotes the derivative with re- spect to ›.

4.6.5 The Local Matrix Equation

Equation 4.42 is suitably described using matrix notation and from now on it will be assumed clear that the integrandsÏÂare functions of › so one may write,

A˙c0

˙c1

B

¸ ˚˙ ˝

˙¯u 1 Õ

ÿ

m=0

ÿ1 n=0

h 2

⁄₁

≠1x²ÏÂmÏÂnd›

¸ ˚˙ ˝

M

= ≠ Ac₀

c₁ B

¸ ˚˙ ˝

¯u

ÿ1 m=0

ÿ1 n=0

2 h

⁄ ₁

≠1x²ÏÂ^Õ_mÏÂ^Õ_nd›,

¸ ˚˙ ˝

K

(4.43)

where the matrices ˙¯u, ¯u, M, and K have been introduced.

¯u is the solution vector where c0 is the value of u at local node ›0and c1is the value of uat local node ›1. ˙¯u is the time derivative of ¯u so that

¯u = Cc0

c₁ D

, ˙¯u =^C˙c0

˙c1

D

. (4.44)

(32)

Matrix M is often referred to as the Mass-matrix (the name comes from early usage of the FEM), and is easily obtained by evaluating the regarded integral for the possible combinations of m and n,

M_m,n=h 2

⁄₁

≠1x²ÏÂ_mÏÂ_nd›. (4.45)

Since m and n each can take 2 values a 2x2 matrix is generated,

M=

CM_0,0 M_0,1 M_1,0 M_1,1 D

, (4.46)

The matrix K is often called the stiffness matrix (also due to early usage of the FEM), the values of K are given by

K_m,n= 2 h

⁄₁

≠1x²ÏÂ^Õ_mÏÂ^Õ_nd›. (4.47)

where m and n may take 2 values each and the possible function values are listed in 4.38. Thus K may be written as a 2x2 matrix

K=

CK_0,0 K_0,1 K_1,0≠ K1,1

D

, (4.48)

Using matrix notation, the FEM diffusion equation may compactly be written as

M ˙¯u = ≠K¯u. (4.49)

In order to find the approximate function u in the entire domain , a global mass matrix, M and a global stiffness matrix K are assembled where the bold font denotes that the matrices are global. The procedure of assembling global matrices is explained in Section 4.8.

26

(33)

4.7. THE REACTION EQUATION

4.7 The reaction equation

In this section the element-matrices describing one-way chemical reactions are derived.

4.7.1 Galerkin’s method Here, the operator L becomes

L(U) = ˆu

dt ≠ ruu, (4.50)

so that the scalar-product of the residual and a test-function v becomes

(R, v) = (ˆu

dt ≠ ruu, v) = (ˆu

dt, v) ≠ (ruu, v). (4.51)

moving the reaction-term to the right hand side and writing the scalar-products in integral form yields

ˆ ˆt

⁄

uvdV = ru

⁄

uvdV, (4.52)

where both integrals have the form of the mass matrix M introduced in the previous section. However, if faced with a reaction of the following kind

On 1

Ó C≠æ U,^r

one sees that the reaction rate depends solely on some constant r and the amount of compound C. If ¯cC(t) and ¯uU(t) denote the solution vectors describing the concentration of compound C and U respectively the equations become in matrix form

C˙cC,0

˙cC,1

D

= ≠rM =

CM_0,0 M_0,1 M_1,0 M_1,1

D Cc_C,0 c_C,1 D

(4.53)

(34)

C˙cU,0

˙cU,1

D

= rM =

CM0,0 M0,1

M_1,0 M_1,1 D Cc_C,0

c_C,1 D

. (4.54)

4.8 Assembling a Global Matrix

In this section local matrices are assembled to create a global matrix, solving for the approximate function u on the entire domain ,

¯u = S WW WW WW U

u₁ u₂ ...

u_N u_N+1

T XX XX XX V

, (4.55)

where the bold font denotes that ¯u is the global solution vector.

As the domain consists of (N + 1) nodes the global matrix is a square matrix with dimension N + 1. Now the element matrix for element n is non-zero on the nodes xn and xn+1 but zero on all other nodes. Note that 2 adjacent elements share a node and therefore overlap in the global matrix. The insertion of the element matrices into a global matrix is illustrated in Figure 4.5.

Once the global matrices are assembled the equation is suitably implemented an in a computational programming environment such as matlab. In this paper the following equations are of special interest,

• The Diffusion Equation: M ˙¯u = ≠K¯u

• The Reaction Equation: M ˙¯u = R¯u

• The Diffusion-Reaction Equation combined: M ˙¯u = (≠K + R)¯u

28

(35)

4.8. ASSEMBLING A GLOBAL MATRIX

Figure 4.5: A global (N + 1) ◊ (N + 1) matrix is assembled by N element matrices.

(36)

(37)

Part II

Model and Method

(38)

(39)

Chapter 5

Model

In this section a mathematical model is established [7], describing what happens to one cell in the in-vitro experiment [2]. First the geometry of the system (consisting of one cell and a finite volume of extra-cellular medium) is established. All local equations acting on the different sub-domains of the system are set up, forming a system of equations.

Further the system is assumed to be spherical symmetrical.

5.1 Geometry

The cell consists of sub-domains listed in Table 2.1 which are presented graphically in Figure 5.1.The entire system is modeled in a perfect spherical fashion, the V79 cells used in the experiment are assumed to have a roughly spherical shape. Further spherical symmetry is assumed on the cell. The membranes are modeled as surfaces, so that they do no posses any volume.

Figure 5.1: The Geometry of a model cell and its surrounding.

(40)

CHAPTER 5. MODEL

5.2 Equations on

i

We seek functions describing the concentration [moles/volume] of the compounds listed in Table ?? throughout the cell. A function describing the concentration for some sub- stance on iis denoted ui.

The substances U and C diffuse and undergo chemical reactions on the aqueous do- mains i, i= 1, 3, 5. Stated mathematically,

ˆu_i(x, t)

ˆt = Ò(Di· Òui(x, t)) + ru,iu_i(x, t), (5.1) For substances B, E and A ,which are produced in chemical reactions but do not diffuse Equation 5.1 simplifies to

ˆui(x, t)

ˆt = ru,i· ui(x, t). (5.2)

5.3 Boundary conditions on

i

As the membranes 1 and 2 are assumed to not possess any volume no diffusion or chemical reactions occur on the membranes, rather substances C and U transmit through the membranes due to being lipophilic. The transmission will occur over the membranes until the concentrations [C] and [U] are at equilibrium in the entire system. Since 1

and 2are infinitely thin the the rate of transmission over a surface is directly proportional to the difference in concentration on the domains adjacent to the regarded surface.

If ni denotes the normal of i, pointing out from the domain according to conven- tion one may describe the transmission over 1as

D1 ˆ

ˆn₁u1= –1(u1≠ u3), (5.3a) D₃ ˆ

ˆn₃u₃= –1(u3≠ u1), (5.3b) (5.3c) where Diis the diffusion constant on i, –1is some proportionality constant. Note also that n1and n3 are oppositely directed, thus

n₁= ≠n3. (5.4)

as illustrated in Figure 5.2. Similarly the transmission over 2is given by

34

(41)

5.4. BOUNDARY CONDITION ON THE OUTER SURFACE

D₃ ˆ

ˆn3u₃= –2(u3≠ u5), (5.5a) D₅ ˆ

ˆn₅u₅= –2(u5≠ u3), (5.5b) where n3and n5are oppositely directed so that

n₁= ≠n5. (5.6)

Substituting Equations 5.4 and 5.6 into Equations 5.3 and 5.5 respectively yields

On 1: D1 ˆ

ˆn₁u₁+ D3 ˆ

ˆn₃u₃= 0, (5.7a)

On 2: D3 ˆ

ˆn₃u₃+ D5 ˆ

ˆn₅u₅= 0. (5.7b)

Equation 5.7b says that what flows out from 1must flow in to 3, and what flows out from 5must flow in to 3, demonstrating that the flux indeed continuous.

Figure 5.2: Normals on respective membranes

5.4 Boundary Condition on the outer surface

In the experimental setup, the volume of the extracellular medium by far exceeds the total volume of all the cells combined, so it is fair to assume there is no exchange of

(42)

CHAPTER 5. MODEL mass between the cells. Thus in our model there is no inflow nor outflow through the constructed surface surrounding the system. In other words, our model consists of a closed system described by the following boundary condition of Neumann-kind,

On : D1 ˆ

ˆn₁u₁= 0. (5.8)

5.5 Lumped parameter approach on the extra-cellular domain,

1

From Equation 5.1 one gathers that partial differential equations (PDEs) 1, describing the concentration of substances C and U,

ˆu1(x, t)

ˆt = Ò(D1Òu1(x, t)) + ru,1· u1(x, t), (5.9)

but these PDEs can be simplified into ordinary differential equations (ODEs) by using a lumped parameters approach. This means that the concentration on the domain is evaluated as if all of the mass was gathered in one place. Using the lumped parameter approach we have assumed the diffusion-rate to be practically instantaneous on the domain. Thus we say 1 is well-stirred and regard the concentration as homogeneous on the whole domain.

The first step to transform the PDEs in Equation 5.9 into ODEs is to write the equations in integral form. Integrating over the volume 1 Equation 5.9 becomes

ˆ ˆt

⁄

1

u₁(x, t)dV = D1

⁄

1 Ò(Òu1(x, t))dV + r1

⁄

1

u₁(x, t)dV, (5.8)

where factors not depending on the spatial location have been brought outside the in- tegrals. Now as u1 denotes the concentration [moles/volume] and the domain is well- stirred, integrating the concentration u1 over the domain-volume yields the number of moles present on 1 denoted Nu,1. Thus the integral equation, Equation 5.8 may be written

ˆ

ˆtN_u,1= D1

⁄

1 Ò(Òu1(x, t))dV + ru,1N_u,1, (5.9)

36

(43)

5.6. LUMPED PARAMETER APPROACH ON THE CELL NUCLEUS, 5

where the remaining integral may be simplified using Gauss’s theorem which can transform volume integral into a closed surface integral. Applying Gauss’s theorem, the remaining volume integral in Equation 5.9 may be rewritten,

D₁

⁄⁄⁄

1Ò · (Òu1)dV =^j

1 (Òu1)dS. (5.10)

With spherical symmetry Equation 5.9 becomes

ˆ

ˆtN_u,1=^j

1

D₁ˆu₁

ˆndS+ ru,1N_u,1. (5.11)

Now using the transmission condition stated in Equation 5.4 a the integrand can be substituted and one obtains

ˆ

ˆtN_u,1=^j

1

–_u,1(u1≠ u3)dS + ru,1N_u,1, (5.12)

evaluating the integral yields

ˆ

ˆtN_u,1= A1–_u,1(u1≠ u3) + r1· Nu,1, (5.13)

where A1denotes the surface-area of 1.

5.6 Lumped parameter approach on the cell nucleus,

5 This domain is treated analogously to 1 as described in Section 5.5 above, therefore the steps taken in this subsection will be described in brief. The equations describing the concentration for Substances C and U are

ˆu₅(x, t)

ˆt = Ò(D5Òu5(x, t)) + ru,5u₅(x, t), (5.14a) (5.14b)

(44)

CHAPTER 5. MODEL The diffusion-rate is assumes to be practically instantaneous so the lumped parameters approach is taken. Integrating Equations 5.14b over the volume 5yields

ˆ ˆt

⁄

5

u₅(x, t)dV = D5

⁄

5 Ò(Òu5(x, t))dV + ru,5

⁄

5

u₅(x, t)dV. (5.15)

If Nu,5 denotes the number of moles on 5 the integral equation, Equation 5.15 becomes

ˆ

ˆtN_u,5= D5

⁄

5 Ò(Òu5(x, t))dV + ru,5N_u,5. (5.16)

Further, applying Gauss’s theorem to the remaining integral in Equation 5.16 above one obtains

ˆ

ˆtNu,5=^j

5

D1ˆu5

ˆndS+ ru,5Nu,5. (5.17)

Using Equation the transmission condition from Equation 5.6 the integrand in Equation 5.17 can be substituted and one obtains

ˆ

ˆtN_u,5=^j

5

–_u,5(u5≠ u3)dS + ru,5N_u,5. (5.18)

evaluating the integral yields

ˆ

ˆtN_u,5= A2–_u,5(u5≠ u3) + ru,5N_u,5, (5.19)

where A2denotes the surface-area of 2.

For substances A and E the diffusion-term in Equation 5.19 disappears and the concentration is described as

38

(45)

5.7. GOVERNING EQUATIONS ON THE CYTOPLASM, 3

ˆ

ˆtN_u,5= ru,5N_u,5. (5.20)

5.7 Governing equations on the cytoplasm,

3

The diffusion-reaction equation governs 3and using boundary conditions obtained from Section 5.3 substances C and U must satisfy the differential equation (D.E) with boundary conditions (B.C)

D.E: ˆ

ˆtu= Ò · (D3Òu3) + r3u₃, (5.21a) B.C on 1: D3ˆu3

ˆn = –1(u3≠ u1), (5.21b) B.C on 2: D3ˆu₃

ˆn = –2(u3≠ u2). (5.21c)

5.8 Evaluating the cytoplasm,

3

The cytoplasm, 3 is complex and highly inhomogeneous therefore it is not suitable to take a lumped-parameters approach when finding functions on this domain. Instead an approximate function u3is obtained using the FEM described in Chapter 4.

Using the weak formulation the diffusion-reaction equation becomes

⁄

3

ˆu₃(x, t) ˆt dV =^⁄

3Ò · (D3Òu3)vdV + r3

⁄

3

u₃vdV, (5.22)

where v is some test function, orthogonal to the solution u3. Partial integration yields the following equation,

⁄

3

ˆu3(x, t)

ˆt u3vdV =^⁄

ˆ 3

D3ˆu3

ˆn₃vdV ≠

⁄

3

D3Òu3Òv3dV + r3

⁄

3

u3vdV, (5.23)

(46)

CHAPTER 5. MODEL

where the first term on the right hand side is a boundary condition. These boundary conditions are given in Equation 5.21, inserting them into Equation 5.23 the governing equation on 3becomes

⁄

3

ˆu₃(x, t)

ˆt vdV = –1

⁄

1

(u3≠u1)vdS+–2

⁄

2

(u3≠u5)vdS≠^⁄

3

D₃Òu3Òv3dV+r3

⁄

3

u₃vdV.

(5.24)

Equation 5.24 is further evaluated using techniques described in Chapter 4.5 where integral equations of this kind is treated.

5.9 Transformation to spherical coordinates

Using spherical coordinates the governing equation on 3becomes

(u3, v) =^⁄

3

u3(x)v(x)dx =^⁄^r^cell

r=rnuc

⁄_2ﬁ

◊=0

⁄_pi

„=0u3(r)v(r)r²sin(„)dr (5.25)

where r dentes radial distance to the center of the cell, „ is the azimuthal angle and

◊ is the polar angle. Now spherical symmetry is assumed on all domains so the scalar product (u3, v) further simplifies to

(u3, v) = 4ﬁ^⁄ ^r^cell

r=rnuc

r²u₃(r)v(r)dr. (5.26)

5.10 Gathering the System of Equations

The total system of equations to solve is now gathered and listed below. A function describing the concentration of substance C on domain i is denoted ciand so on for all

40

(47)

5.10. GATHERING THE SYSTEM OF EQUATIONS compounds.

V₅d

dtc₅= ≠–2A₂(c5≠ c3(ﬂnuc) + V5r_C,5c₅ (5.27a) V₅d

dtu₅= ≠–2A₂(u5≠ u3(ﬂnuc) + V5r_U,5c₅ (5.27b) V₅d

dta₅= V5r_A,5c₅ (5.27c)

V₅d

dte₅= V5r_E,5c₅ (5.27d)

V1d

dtc1= ≠–1A1(c1≠ c3(ﬂcell) + V1rC,1c1 (5.27e) V₁d

dtu₁= ≠–1A₁(u1≠ u3(ﬂcell) + V1r_U,5c₁ (5.27f) d

dt(u3, v) = ≠–2A₂(c5≠ c3(flnuc)v(flnuc) ≠ –1A₁(c1≠ c3(flcell)v(flcell) ≠ (D3Òc3,Òv) + (rC,3c₃, v) (5.27g) d

dt(u3, v) = ≠–2A2(u5≠ u3(flnuc)v(flnuc) ≠ –1A1(u1≠ u3(flcell)v(flcell) ≠ (D3Òu3,Òv) + (rU,3c3, v) (5.27h)

d

dt(b3, v) = (rB,3c₃, v) (5.27i)

d

dt(e3, v) = (rE,3c₃, v) (5.27j)

(48)

(49)

Chapter 6

Method

This section provides a method for solving the system of equations, Equations 5.27.

6.1 Discretization of the entire domain

The functions on the cytoplasm, 3, are evaluated using the FEM, thus the domain 3is discretized and contains N elements and therefore N + 1 nodes as explained in Chapter 4. The discretization is now expanded to include the entire cell. Note that the functions on 5 are not dependent on spatial location so one node is sufficient to represent the function-value on the entire domain 5, the same is true for function-values on 1. Thus if x0is a node on 5and xN+2is a node on 1 the discretized domain representing the entire cell comprises the nodes x0, x₁, ...x_N+1, x_N+2 as illustrated in Figure 6.1. Note especially that x1 lies on the nucleus membrane where r = ﬂnucand xN+1 lies on the cell membrane where r = ﬂcell.

Figure 6.1: The domain is discretized and consists of N + 2 elements if 3consists of N elements. The node x0lies in 5, node x1lies on 2, node xN+1lies on 1and node x_N+2 lies in 1.

(50)

CHAPTER 6. METHOD

6.2 A matrix-equation for the entire system of equations

To solve for all unknowns, that is all function-values the different functions (c, u, b, e, a) take on each node, a solution vector ¯w is formed,

¯ w =

Q cc cc cc cc cc cc cc cc cc cc cc cc cc cc ca

c₀ u₀ a₀ e₀ c1

u₁ b₁ e₁ ...

cN+1

u_N+1 b_N+1 e_N₊₁ c_N+2 uN+2

R dd dd dd dd dd dd dd dd dd dd dd dd dd dd db

, (6.1)

where cn denotes the funtion-value c takes on node xn and so on. To solve the system of equations the FEM matrix-form is used so that one has an equation of the form

M_full˙¯w = (≠Kfull+ R_full)¯w, (6.2)

where M is the mass-matrix and K is the stiffness-matrix as described in Section 4.6.5.

The subscripts full denote that the matrices represent all nodes on .

Disregarding the boundary conditions for now, each full matrix is a block-matrix where the domains 1, 3, 5 are represented by one block each. Figure ?? illustrates the block-matrix form of some equation ˙¯w = A¯w, for some matrix A. The explicit expres- sions for the matrix-blocks are given in Subsections 6.2.1, 6.5, 6.2.3.

44

(51)

6.2. A MATRIX-EQUATION FOR THE ENTIRE SYSTEM OF EQUATIONS

Figure 6.2: A block-matrix representation of the equation ˙¯w = A¯w, where ¯w is a (N + 3) ◊ 1 column vector and A is a (N + 3) ◊ (N + 3) square-matrix.

6.2.1 Matrices on 5

The matrices on are

¯ w ₅=

Q cc ca c₀ u₀ a₀ e₀

R dd db,M ₅ =

Q cc ca

1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 R dd db,K ₅ =

Q cc ca

D_C,3 0 0 0 0 D_U,3 0 0

0 0 0 0

R dd db,R ₅ =

Q cc ca

r_C,5 0 0 0 r_U,5 1 0 0 r_A,5 0 0 0 r_E,5 0 0 0 R dd db,

(6.3)

6.2.2 Matrices on 1

The matrices on 1are

¯ w ₁=

A c_N₊₂ u₀N+ 2

B

,M ₁= A1 0

0 1 B

,K ₁ =

AD_C,1 0 0 D_U,1

B ,R ₅=

Ar_C,1 0 r_U,1 0 B

. (6.4)

A numerical approach to model the dynamics in a mammalian cell exposed to a carcinogenic compound