1002/1002.4976.md

Analytical And Numerical Approximation of Effective Diffusivities in The Cytoplasm of Biological Cells

Michael Hanke and Marry-Chriz Cabauatan-Villanueva

Royal Institute of Technology

School of Computer Science and Communication

June 11, 2007

Abstract

The simulation of the metabolism in mammalian cells becomes a severe problem if spatial distributions must be taken into account. Especially the cytoplasm has a very complex geometric structure which cannot be handled by standard discretization techniques. In the present paper we propose a homogenization technique for computing effective diffusion constants. This is accomplished by using a two-step strategy. The first step consists of an analytic homogenization from the smallest to an intermediate scale. The homogenization error is estimated by comparing the analytic diffusion constant with a numerical estimate obtained by using real cell geometries. The second step consists of a random homogenization. Since no analytical solution is known to this homogenization problem, a numerical approximation algorithm is proposed. Although rather expensive this algorithm provides a reasonable estimate of the homogenized diffusion constant.

Contents

1	Introduction	2
2	The Mathematical Model	4
2.1	The Governing Equations . . . . .	4
2.2	The Model Problem . . . . .	5
2.3	Going From The Smallest to The Medium Scale: Homogenization of a Periodic Structure . . . . .	6
2.4	From The Medium to The Large Scale: Stochastic Homogenization . . . . .	8
3	Numerical Determination of Effective Diffusivities	10

4	Theoretical And Experimental Diffusivities For Layered Structures	11
4.1	The Experimental Set-up . . . . .	12
4.2	Results . . . . .	14
5	Experimental Effective Diffusivities in Random Media	16
5.1	The Experimental Set-up . . . . .	16
5.2	Results in 2D . . . . .	17
5.3	Results in 3D . . . . .	18
6	Conclusions	18

1 Introduction

When mammalian cells are exposed to foreign and potentially harmful compounds a series of events takes place. Following uptake the substance is distributed in different intracellular compartments by diffusion, absorption and desorption. The majority of the compound is either dissolved in the aqueous phase, the cytoplasm, or in the lipophilic phase, the membranes. Parallel to diffusion and absorption/desorption bioactivation/biotransformation by different soluble and membrane bound enzymes takes place. The purpose of biotransformation is to render the substance suitable for excretion.

A human cell consists schematically of an outer cellular membrane, a cytoplasm containing a large number of organelles (mitochondria, endoplasmatic reticulum etc.), a nuclear membrane and finally the cellular nucleus containing DNA. Figure 1 shows a sketch of a cell while Figure 2 shows a microphotograph of a nucleus with part of the surrounding cytoplasm. The organelle membranes create a complex and dense system of membranes or subdomains throughout the cytoplasm. The mathematical description of the biotransformation leads to a system of reaction-diffusion equations in a complex geometrical domain, dominated by thin membranous structures with similar physical and chemical properties. If these structures are treated as separate subdomains, any model becomes computationally very expensive. Moreover, due to the natural variations in the cell structures, every individual cell needs its own mathematical model.

In order to make the system numerically treatable while at the same time retaining the essential features of the metabolism under consideration, in [4] a way of homogenizing the cytoplasm has been developed, aiming at a manageable system of reaction-diffusion equations for the various species. In the present paper, we report about numerical experiments which justify some of the strategies in the cited paper. The general modelling assumptions are summarized below. We will use them also in the present report.

Modelling assumptions:

• • On a small scale in space, the volume between the outer cellular membrane and the nucleus membrane consists of layered structures cytoplasm/membranes.
• • In the large scale, this volume contains an unordered set of the small-scale substructures which are uniformly distributed over the volume.Figure 1: Schematic picture of a cell. Picture copyright U.S. National Cancer Institute's Surveillance, Epidemiology and End Results Program, http://training.seer.cancer.gov/module\_anatomy/unit2\_1\_cell\_functions\_1.html

Figure 2: Ultrastructure of the cell, nucleus and cytoplasm. Picture copyright Histology Learning System, Boston University, http://www.bu.edu/histology/m/index.htm- • The physical and chemical properties of the cytoplasm and of the membranes are uniform.

• • We adopt the continuum hypothesis, i.e., we assume that the set of molecules in the cell can be modelled by considering a continuous representation (a concentration).
• • The processes of absorption and desorption of the individual species into or out of the membrane is much faster than the diffusion and reaction processes. In this case, the relations of the concentrations of a species near a membrane/cytoplasm boundary can be conveniently described with the help of a partition coefficient.

In Section 2, we introduce the mathematical model. The next section is devoted to a description of the general experimental set-up which is used to compute effective diffusion coefficients numerically. For the solution of the arising boundary-value problems for partial differential equations we used the Comsol Multiphysics¹ [2] environment.

In [?], the diffusion coefficients in the membrane structures have been homogenized by assuming the membranes and the aqueous volumes to be ideal infinite plane layers. This allows for an analytic computation of the effective diffusivity. In Section 4, we will compare effective diffusivities obtained this way against numerically determined effective diffusivities by using computational domains which have been discretized from microphotographs of cell membranes.

The result of the first homogenization step leads to anisotropic diffusion tensors valid locally. Invoking the assumption about the random distribution of the orientation of the membranes, the next step consists of a stochastic homogenization. In contrast to the one- and two-dimensional case, no analytical solution in the general three-dimensional case is known. We will compute the effective diffusion coefficient numerically by Monte Carlo techniques in Section 5.

2 The Mathematical Model

2.1 The Governing Equations

We intend to derive a homogenized model of the reaction and diffusion processes inside the cytoplasm. For that purpose, let $G$ denote the volume between the outer cell membrane and the nucleus membrane (excluding the membranes themselves). This volume is split into two disjoint parts $G_l$ and $G_w$ which denote the lipophilic part and the aqueous part, respectively, of the cell. Note that these subdomains are not necessarily connected. Assume that we are interested in the concentrations $c_1, \dots, c_n$ of $n$ species inside of $G$ . For the $k$ -th species, it holds

$$\frac{\mathrm{\partial }}{\mathrm{\partial }t}{c}_k=\mathrm{\nabla }\cdot (d_k(x)c_k)+R_k(c_1,\dots ,c_n,x),x\in G,k=1,\dots ,n\mathrm{.}(1)\backslash frac\{\backslash partial\}\{\backslash partial\; t\}\; c\_k\; =\; \backslash nabla\; \backslash cdot\; (d\_k(x)\; c\_k)\; +\; R\_k(c\_1,\; \backslash dots,\; c\_n,\; x),\; \backslash quad\; x\; \backslash in\; G,\; \backslash quad\; k\; =\; 1,\; \backslash dots,\; n.\; \backslash quad\; (1)$$

Here, $d_k$ denotes the diffusion tensor of the $k$ -th species which is assumed to be constant in both $G_l$ and $G_w$ . $R_k$ denotes the reaction term. It varies strongly with $x$ . In the lipophilic part, $R_k \equiv 0$ because no reactions are taking place there. The concentrations of some of the species can be assumed to be constant over time. As a consequence, many of the reaction terms will be linear.

---

¹Comsol Multiphysics is a registered trade mark of Comsol AB, Stockholm, Sweden.The partition coefficient, $K_{p,k}$ is the equilibrium ratio of the concentrations of species $k$ between the aqueous compartment and its adjacent lipid compartment. This gives rise to boundary conditions,

$$c_{k,w}=K_{p,k}{c}_{k,l}x\in G_w\cap G_l,(2)c\_\{k,w\}\; =\; K\_\{p,k\}\; c\_\{k,l\}\; \backslash quad\; x\; \backslash in\; G\_w\; \backslash cap\; G\_l,\; \backslash quad\; (2)$$

on the inner boundaries where $c_{k,w}$ and $c_{k,l}$ denote the concentrations in the aqueous and lipid parts, respectively.

The system (1) with inner boundary conditions (2) will be supplemented by (outer) boundary conditions and initial conditions. For the purposes of this paper the precise structure of these conditions is not important.

2.2 The Model Problem

Let $G \subset \mathbb{R}^3$ be a bounded domain which will be splitted into two (not necessarily connected) subdomains $G_1$ and $G_2$ such that

$$G_1\cap G_2=\mathrm{\varnothing },\bar{G}={\bar{G}}_1\cup {\bar{G}}_2\mathrm{.}G\_1\; \backslash cap\; G\_2\; =\; \backslash emptyset,\; \backslash quad\; \backslash bar\{G\}\; =\; \backslash bar\{G\}\_1\; \backslash cup\; \backslash bar\{G\}\_2.$$

The interior boundary will be denoted by $\Gamma$ ,

$$\mathrm{\Gamma }={\bar{G}}_1\cap {\bar{G}}_2\mathrm{.}\backslash Gamma\; =\; \backslash bar\{G\}\_1\; \backslash cap\; \backslash bar\{G\}\_2.$$

Consider the equations

$$\frac{\mathrm{\partial }}{\mathrm{\partial }t}{v}_i-\mathrm{\nabla }\cdot (d_i(x)\mathrm{\nabla }{v}_i)+r_i(x)v_i=f_i(x),x\in G_i,i=1,2.(3)\backslash frac\{\backslash partial\}\{\backslash partial\; t\}\; v\_i\; -\; \backslash nabla\; \backslash cdot\; (d\_i(x)\; \backslash nabla\; v\_i)\; +\; r\_i(x)\; v\_i\; =\; f\_i(x),\; \backslash quad\; x\; \backslash in\; G\_i,\; \backslash quad\; i\; =\; 1,\; 2.\; \backslash quad\; (3)$$

Assume additionally boundary conditions on $\partial G$ and initial conditions on $G$ be given.

On the inner boundary $\Gamma$ the flux must be continuous. Let $\mathbf{n}_i$ denote the outer normal at the boundary of $G_i$ . The continuity conditions reads now

$$d_1\frac{\mathrm{\partial }{v}_1}{\mathrm{\partial }{\mathbf{n}}_1}+d_2\frac{\mathrm{\partial }{v}_2}{\mathrm{\partial }{\mathbf{n}}_2}=0,x\in \mathrm{\Gamma }\mathrm{.}(4)d\_1\; \backslash frac\{\backslash partial\; v\_1\}\{\backslash partial\; \backslash mathbf\{n\}\_1\}\; +\; d\_2\; \backslash frac\{\backslash partial\; v\_2\}\{\backslash partial\; \backslash mathbf\{n\}\_2\}\; =\; 0,\; \backslash quad\; x\; \backslash in\; \backslash Gamma.\; \backslash quad\; (4)$$

The presence of a partion coefficient between the two phases gives rise to the boundary condition

$$v_1=K_p{v}_2,x\in \mathrm{\Gamma }\mathrm{.}(5)v\_1\; =\; K\_p\; v\_2,\; \backslash quad\; x\; \backslash in\; \backslash Gamma.\; \backslash quad\; (5)$$

This problem can be reduced to a problem in a more standard form by introducing

For this new function $u$ , the inner boundary conditions become

This motivates the definitions

With these definitions, the problem (4), (5) becomes equivalent to

$$\sigma \frac{\mathrm{\partial }}{\mathrm{\partial }t}u-\mathrm{\nabla }\cdot (d\mathrm{\nabla }u)+ru=f,x\in G(7)\backslash sigma\; \backslash frac\{\backslash partial\}\{\backslash partial\; t\}\; u\; -\; \backslash nabla\; \backslash cdot\; (d\; \backslash nabla\; u)\; +\; ru\; =\; f,\; \backslash quad\; x\; \backslash in\; G\; \backslash quad\; (7)$$

subject to correspondingly modified initial and boundary conditions.

For later use, let

$$p_1=\frac{\mathrm{\mid }{G}_1\mathrm{\mid }}{\mathrm{\mid }G\mathrm{\mid }},p_2=1-p_1=\frac{\mathrm{\mid }{G}_2\mathrm{\mid }}{\mathrm{\mid }G\mathrm{\mid }}\mathrm{.}(8)p\_1\; =\; \backslash frac\{|G\_1|\}\{|G|\},\; \backslash quad\; p\_2\; =\; 1\; -\; p\_1\; =\; \backslash frac\{|G\_2|\}\{|G|\}.\; \backslash quad\; (8)$$

2.3 Going From The Smallest to The Medium Scale: Homogenization of a Periodic Structure

The homogenization procedure for an equation of the type (7) is proved in [6]. We cite the basic facts. Consider the following problem,

Here, the operator $A^\varepsilon$ is given by

$$A^{\epsilon }u:=-\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(d_{ij}^{\epsilon }\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}u\right)+r^{\epsilon }u\mathrm{.}A^\backslash varepsilon\; u\; :=\; -\backslash frac\{\backslash partial\}\{\backslash partial\; x\_i\}\; \backslash left(\; d\_\{ij\}^\backslash varepsilon\; \backslash frac\{\backslash partial\}\{\backslash partial\; x\_j\}\; u\; \backslash right)\; +\; r^\backslash varepsilon\; u.$$

For convenience, we use the Einstein summation convention: If an index appears twice in a multiplicative expression, this expression is understood to implicitly represent the sum over this expression where the index varies between 1 and 3 (the dimension of $G$ ). Moreover, we assume the following construction of the coefficients:

$${\sigma }^{\epsilon }(x)=\sigma (x\mathrm{/}\epsilon ),r^{\epsilon }(x)=r(x\mathrm{/}\epsilon ),d_{ij}^{\epsilon }(x)=d_{ij}(x\mathrm{/}\epsilon ),i,j=1,2,3.\backslash sigma^\backslash varepsilon(x)\; =\; \backslash sigma(x/\backslash varepsilon),\; \backslash quad\; r^\backslash varepsilon(x)\; =\; r(x/\backslash varepsilon),\; \backslash quad\; d\_\{ij\}^\backslash varepsilon(x)\; =\; d\_\{ij\}(x/\backslash varepsilon),\; \backslash quad\; i,\; j\; =\; 1,\; 2,\; 3.$$

The functions $\sigma^\varepsilon$ and $r^\varepsilon$ are assumed to belong to $L^\infty(G)$ , and

$$\sigma \ge {\sigma }_0>0,r(x)\ge 0\text{ a.e. in }G\backslash sigma\; \backslash geq\; \backslash sigma\_0\; >\; 0,\; \backslash quad\; r(x)\; \backslash geq\; 0\; \backslash text\{\; a.e.\; in\; \}\; G$$

for some $\sigma_0 \in \mathbb{R}$ . The functions $d_{ij}^\varepsilon$ are assumed to be measurable and to satisfy the conditions $d_{ij}^\varepsilon = d_{ji}^\varepsilon$ and

### 2.3 Going From The Smallest to The Medium Scale: Homogenization of a Periodic Structure⁷

Finally, assume $u_0 \in L^2(G)$ and $f^\varepsilon \in L^2(0, T; L^2(G))$ .

In order to find the homogenized equation, assume

$$f^{\epsilon }⟶f\text{ weakly in }{L}^2(0,T;L^2(G))\mathrm{.}f^\backslash varepsilon\; \backslash longrightarrow\; f\; \backslash text\{\; weakly\; in\; \}\; L^2(0,\; T;\; L^2(G)).$$

Let $Y$ be an axis parallel hexahedron in $\mathbb{R}^3$ , that is,

$$Y=(a_i,b_i)\mathrm{.}Y\; =\; \backslash underset\{i=1\}\{\backslash overset\{3\}\{\backslash times\}\}\; (a\_i,\; b\_i).$$

For a $Y$ -periodic function $f$ , the mean value is given by

$$⟨f⟩:=\frac{1}{\mathrm{\mid }Y\mathrm{\mid }}{\int }_{Y}f(y)dy\mathrm{.}\backslash langle\; f\; \backslash rangle\; :=\; \backslash frac\{1\}\{|Y|\}\; \backslash int\_Y\; f(y)\; dy.$$

Assume now that $a_{ij}$ , $\sigma$ , and $r$ are $Y$ -periodic. Then it is possible to consider the problem

where the operator $A$ is given by

and $\varphi_j$ is the $Y$ -periodic solution of the following local elliptic problem:

Here, $W(Y) = {\phi \in H^1(Y) \mid \phi \text{ is } Y\text{-periodic and } \langle \phi \rangle = 0}$ .

Theorem 2.1. Under the conditions stated above, (9) and (10) have unique solutions $u^\varepsilon \in L^2(0, T; H_0^1(G))$ and $u \in L^2(0, T; H_0^1(G))$ , respectively, and it holds

$$u^{\epsilon }⟶u\text{ in }{L}^2(0,T;H_0^1(G))\text{ weakly as }\epsilon \to 0.u^\backslash varepsilon\; \backslash longrightarrow\; u\; \backslash text\{\; in\; \}\; L^2(0,\; T;\; H\_0^1(G))\; \backslash text\{\; weakly\; as\; \}\; \backslash varepsilon\; \backslash rightarrow\; 0.$$

This theorem is proved in [6, p. 56].²

In our model problem (7) the cell problems can be simplified considerably. We will assume that, in the smallest scale, aqueous and lipid compartments are perfectly layered. It turns out

---

²In the reference, the proof is given for a problem without reaction term. But the proof can easily be generalized.that in this case the computation for the effective diffusion coefficients leads to a transmission problem which can be solved analytically.

Consider the cell problem (11). The material consists of perfect layers $\omega^+$ and $\omega^-$ with thicknesses $a^+$ and $a^-$ , respectively. The diffusion coefficients in these layers are $d^+$ and $d^-$ , respectively. According to our assumptions,

$$\frac{a^{+}}{a^{-}}=\frac{p_1}{p_2}\mathrm{.}\backslash frac\{a^+\}\{a^-\}\; =\; \backslash frac\{p\_1\}\{p\_2\}.$$

Let us choose the following coordinate: $x_2$ is normal to the interface between $\omega^+$ and $\omega^-$ while $x_1$ and $x_3$ span this interface. The diffusion coefficient $d$ is given by

The cell problem is posed on

$$Y=(0,l_1)\times (-a^{-},a^{+})\times (0,l_3)\mathrm{.}Y\; =\; (0,\; l\_1)\; \backslash times\; (-a^-,\; a^+)\; \backslash times\; (0,\; l\_3).$$

Then the homogenized diffusivities become (see [7])

$$d_{\text{eff},11}=(a^{+}{d}_1^{+}+a^{-}{d}_1^{-})\mathrm{/}(a^{+}+a^{-}),(12)d\_\{\backslash text\{eff\},11\}\; =\; (a^+\; d\_1^+\; +\; a^-\; d\_1^-)\; /\; (a^+\; +\; a^-),\; \backslash quad\; (12)$$

$$d_{\text{eff},22}=(a^{+}+a^{-})\mathrm{/}(a^{+}\mathrm{/}{d}_2^{+}+a^{-}\mathrm{/}{d}_2^{-}),(13)d\_\{\backslash text\{eff\},22\}\; =\; (a^+\; +\; a^-)\; /\; (a^+\; /\; d\_2^+\; +\; a^-\; /\; d\_2^-),\; \backslash quad\; (13)$$

$$d_{\text{eff},33}=(a^{+}{d}_3^{+}+a^{-}{d}_3^{-})\mathrm{/}(a^{+}+a^{-}),(14)d\_\{\backslash text\{eff\},33\}\; =\; (a^+\; d\_3^+\; +\; a^-\; d\_3^-)\; /\; (a^+\; +\; a^-),\; \backslash quad\; (14)$$

$$d_{ij}=0\text{ if }i\mathrm{\ne }j\mathrm{.}(15)d\_\{ij\}\; =\; 0\; \backslash text\{\; if\; \}\; i\; \backslash neq\; j.\; \backslash quad\; (15)$$

Note that $d_{\text{eff},11}$ and $d_{\text{eff},33}$ are the arithmetic means while $d_{\text{eff},22}$ is the harmonic mean of both diffusivities $a_i^+$ and $a_i^-$ .

2.4 From The Medium to The Large Scale: Stochastic Homogenization

In global coordinates, we cannot assume that the coordinate system is oriented in the way that we used above. Consider two Cartesian coordinate systems $(x_1, x_2, x_3)$ and $(z_1, z_2, z_3)$ . Assume that a given point $x$ has the representation $z = Tx$ with respect to the $z$ -coordinates. Note that $T$ is an orthogonal matrix in that case. Denote the matrix of diffusion coefficients with respect to the $x$ -coordinates by $Q^*$ and that with respect to the $z$ -coordinates by $Q$ . Then a short calculation yields

$$Q^{\ast }=TQ{T}^{-1}=TQ{T}^T\mathrm{.}Q^*\; =\; TQT^\{-1\}\; =\; TQT^T.$$

This is the point to invoke the next critical assumption: We assume that the volume is tightly packed with substructures of the type considered before, namely layered materials. The key assumption is that all orientations are equally probable. This leads to a stochastic description ofthe diffusion coefficients. We need a homogenization of operators with random coefficients. A theory for that is provided in [5].

Note that the mean values $\langle \sigma \rangle$ and $\langle r \rangle$ in (10) are independent of the orientation of the layers. Therefore, it suffices to consider the stationary diffusion problem

$$A^{\epsilon }{u}^{\epsilon }=f,u^{\epsilon }\in H_0^1(G)(16)A^\backslash varepsilon\; u^\backslash varepsilon\; =\; f,\; \backslash quad\; u^\backslash varepsilon\; \backslash in\; H\_0^1(G)\; \backslash quad\; (16)$$

with $G \subset \mathbb{R}^m$ and

$$A^{\epsilon }u=-\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left(d_{ij}^{\epsilon }\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}u\right)A^\backslash varepsilon\; u\; =\; -\backslash frac\{\backslash partial\}\{\backslash partial\; x\_i\}\; \backslash left(\; d\_\{ij\}^\backslash varepsilon\; \backslash frac\{\backslash partial\}\{\backslash partial\; x\_j\}\; u\; \backslash right)$$

which is the counterpart of (9). Assume as before that

$$d_{ij}^{\epsilon }(x)=d_{ij}(x\mathrm{/}\epsilon ),i,j=1,\dots ,m\mathrm{.}d\_\{ij\}^\backslash varepsilon(x)\; =\; d\_\{ij\}(x/\backslash varepsilon),\; \backslash quad\; i,\; j\; =\; 1,\; \backslash dots,\; m.$$

The randomness of the orientation is modelled by assuming that the matrix $A(y) = (d_{ij}(y))$ is statistically stationary with respect to the spatial variable $y \in \mathbb{R}^m$ , or equivalently, that $A(y)$ is a typical realization of a stationary random field.

Let $(\Omega, \mathcal{F}, P)$ be a probability space with $\sigma$ -algebra $\mathcal{F}$ and probability measure $P$ . Let for each $x \in \mathbb{R}^m$ a random variable $\xi(x)$ over $(\Omega, \mathcal{F}, P)$ be given. The random field $\xi$ is stationary if it can be represented in the form

$$\xi (x,\omega )=a(T(x)\omega )\backslash xi(x,\; \backslash omega)\; =\; a(T(x)\backslash omega)$$

where $a(\cdot)$ is a fixed random variable, $T = T(x) : \Omega \rightarrow \Omega$ is a measurable transformation which preserves the measure $P$ on $(\Omega, \mathcal{F})$ . Therefore, for the definition of the coefficients $d_{ij}$ in (16) it is sufficient to consider a matrix $(d_{ij})$ of random variables $d_{ij} : \Omega \rightarrow \mathbb{R}$ . Realizations of coefficients can then be obtained by setting

$$d_{ij}(x,\omega )=d_{ij}(T(x)\omega )\mathrm{.}d\_\{ij\}(x,\; \backslash omega)\; =\; d\_\{ij\}(T(x)\backslash omega).$$

Assume in the following that $d_{ij} \in L^\infty(\Omega)$ and

$$\alpha \mathrm{\mid }\xi {\mathrm{\mid }}^2\le d_{ij}(\omega ){\xi }_i{\xi }_j,\xi \in {\mathbb{R}}^m\backslash alpha\; |\backslash xi|^2\; \backslash leq\; d\_\{ij\}(\backslash omega)\; \backslash xi\_i\; \backslash xi\_j,\; \backslash quad\; \backslash xi\; \backslash in\; \backslash mathbb\{R\}^m$$

for almost all $\omega \in \Omega$ with $\alpha > 0$ independent of $\xi$ and $\omega$ .

A (deterministic) matrix $d(y)$ is said to admit a homogenization if there exists a constant elliptic matrix $d_{\text{eff}}$ such that for any $f \in H^{-1}(G)$ the solutions $u^\varepsilon$ of the Dirichlet problem (16) it holds

where $u$ is the solution of the Dirichlet problem

$$-\mathrm{\nabla }\cdot (d_{\text{eff}}\mathrm{\nabla }u)=f,u\in H_0^1(G)\mathrm{.}-\backslash nabla\; \backslash cdot\; (d\_\{\backslash text\{eff\}\}\; \backslash nabla\; u)\; =\; f,\; \backslash quad\; u\; \backslash in\; H\_0^1(G).$$

This definition corresponds to the stationary version of Theorem 2.1. The following theorem holds true [5, p. 230]:Theorem 2.2. Assume additionally that the family of mappings $T(x) : \Omega \rightarrow \Omega$ , $x \in \mathbb{R}^m$ , forms an ergodic $m$ -dimensional dynamical system. Then for almost all $\omega \in \Omega$ , the matrix with coefficients $d_{ij}(x) = d_{ij}(T(x)\omega)$ admits homogenization, and the homogenized matrix $d_{\text{eff}}$ is independent of $\omega$ .

Unfortunately, an analytical representation of $d_{\text{eff}}$ is only possible in exceptional cases. We are interested in diffusion coefficients having a representation

$$d(x)=T(x,\omega )QT(x,\omega {)}^{-1}d(x)\; =\; T(x,\; \backslash omega)QT(x,\; \backslash omega)^\{-1\}$$

where $T(x) \in SO(m)$ is uniformly distributed in $SO(m)$ and $Q$ is a fixed diffusion tensor. In the two-dimensional case, an analytic solution is provided in [5, p. 235]. For $m = 2$ , $d_{\text{eff}}$ is simply a scalar equal to the geometric mean of the eigenvalues of $Q$ ,

$$d_{\text{eff}}=\sqrt{\det (Q)}\mathrm{.}d\_\{\backslash text\{eff\}\}\; =\; \backslash sqrt\{\backslash det(Q)\}.$$

Here $\det(Q)$ denotes the determinant of $Q$ .

There is no analytical solution known for the case $m = 3$ .

For later use in the experimental estimation of the effective diffusivity, the following observation is important: According to our assumptions on $d$ , the estimate

$$\mathrm{\parallel }(A^{\epsilon }{)}^{-1}\mathrm{\parallel }\le {\alpha }^{-1}\backslash |(A^\backslash varepsilon)^\{-1\}\backslash |\; \backslash leq\; \backslash alpha^\{-1\}$$

holds true such that, for any $f$ and $l$ in $H^{-1}(G)$ ,

$$\mathrm{\mid }⟨l,(A^{\epsilon }{)}^{-1}f⟩\mathrm{\mid }\le \mathrm{\parallel }l{\mathrm{\parallel }}_{H^{-1}(G)}{\alpha }^{-1}\mathrm{\parallel }f{\mathrm{\parallel }}_{H^{-1}(G)}|\backslash langle\; l,\; (A^\backslash varepsilon)^\{-1\}\; f\; \backslash rangle|\; \backslash leq\; \backslash |l\backslash |\_\{H^\{-1\}(G)\}\; \backslash alpha^\{-1\}\; \backslash |f\backslash |\_\{H^\{-1\}(G)\}$$

independently of $\omega \in \Omega$ . Consequently, for the expectation values it holds

$$\mathbb{E}⟨l,(A^{\epsilon }{)}^{-1}f⟩⟶\mathbb{E}⟨l,A^{-1}f⟩(17)\backslash mathbb\{E\}\backslash langle\; l,\; (A^\backslash varepsilon)^\{-1\}\; f\; \backslash rangle\; \backslash longrightarrow\; \backslash mathbb\{E\}\backslash langle\; l,\; A^\{-1\}\; f\; \backslash rangle\; \backslash quad\; (17)$$

by the dominated convergence theorem.

3 Numerical Determination of Effective Diffusivities

Under the assumption that an effective diffusivity for a given problem exists, the corresponding diffusion constants can be determined experimentally. For that, let $D \subset G$ be a subdomain which is in size comparable to $G$ such that the small scale structure is considerably smaller than the size of $D$ . Assume that we want to determine the (scalar) diffusion constant for the diffusion process in $x$ -direction. In that case it is convenient to use a cylindrical domain

$$D=(0,L)\times \omega D\; =\; (0,\; L)\; \backslash times\; \backslash omega$$

with $\omega \subset \mathbb{R}^2$ being some bounded domain. On $D$ consider the stationary diffusion equation

$$-\mathrm{\nabla }\cdot (d(x)\mathrm{\nabla }u)=0,x\in D\mathrm{.}-\backslash nabla\; \backslash cdot\; (d(x)\backslash nabla\; u)\; =\; 0,\; \backslash quad\; x\; \backslash in\; D.$$

The boundary conditions are selected as follows:- • On the boundary $\Gamma_0 = {0} \times \omega$ , a fixed Dirichlet condition is given,

$$u_{{\mathrm{\Gamma }}_0}=c_0\mathrm{.}u\_\{\backslash Gamma\_0\}\; =\; c\_0.$$

• • On the boundary $\Gamma_1 = {1} \times \omega$ , a free diffusion into the surrounding medium is assumed,

$$-n\cdot (d(x)\mathrm{\nabla }u{)}_{{\mathrm{\Gamma }}_1}=M(u_{{\mathrm{\Gamma }}_1}-c_1)\mathrm{.}-n\; \backslash cdot\; (d(x)\; \backslash nabla\; u)\_\{\backslash Gamma\_1\}\; =\; M(u\_\{\backslash Gamma\_1\}\; -\; c\_1).$$

Here, $M$ is the mass transfer coefficient and $c_1$ is the concentration in the bulk solution outside of $D$ .

• • All other boundaries $\Gamma_2 = \partial G \setminus (\Gamma_0 \cup \Gamma_1)$ are isolated,

$$-n\cdot (d(x)\mathrm{\nabla }u{)}_{{\mathrm{\Gamma }}_2}=0.-n\; \backslash cdot\; (d(x)\; \backslash nabla\; u)\_\{\backslash Gamma\_2\}\; =\; 0.$$

If $d(x)$ would be a constant $d_{\text{eff}}$ , it would hold

By $|\Gamma_1|$ we denote the Lebesgue measure of $\Gamma_1$ . If $d(x)$ is varying, these equations can be used as an estimation of the effective diffusivity $d_{\text{eff}}$ . In case of an anisotropic effective diffusivity, the above construction leads to an estimate of the effective diffusivity in $x$ -direction, i.e., $d_{\text{eff},11}$ .

In the one-dimensional case $m = 1$ , however, an analytic solution is possible. A simple calculation gives

$$d_{\text{eff}}={(\frac{1}{L}{\int }_0^{L}d(x{)}^{-1}dx)}^{-1},d\_\{\backslash text\{eff\}\}\; =\; \backslash left(\; \backslash frac\{1\}\{L\}\; \backslash int\_0^L\; d(x)^\{-1\}\; dx\; \backslash right)^\{-1\},$$

which amounts to the harmonic mean.

4 Theoretical And Experimental Diffusivities For Layered Structures

The homogenization of layered structures in Section 2.3 made use of the assumption that we have ideal planes of different materials with different diffusion tensors. In a real biological cell, this assumption is only approximately fulfilled in small subdomains. Besides the effect of not having the parameter $\varepsilon$ close to zero an additional error is introduced this way. The aim of the present section is to obtain some experimental estimates of how large the error will be. We will start with a real photograph of some cell organelles and extract the geometrical structure of the lipophilic and aqueous layers. Then the diffusivity is estimated using the strategy of Section 3. This diffusion constant will be compared to the theoretical homogenized value.Figure 3: Detail of a rat cell showing the Golgi-apparatus. The box indicates the area used as a reference domain. Copyright Dr. H. Jastrow

4.1 The Experimental Set-up

The experiments of this section are based on the micro-photograph shown in Figure 3. The part enclosed by a box in that figure has been extracted and amplified in contrast. This way, the membrane structure in Figure 4 has been obtained. Note that only the black lines represent membranes. The geometry of this structure has been too complex for the software used in the numerical experiments. The number of degrees of freedoms obtained after discretization has become too large. Therefore, we extracted again a part of this geometry in order to make the problem tractable with the available software. Note that the diffusion in this problem is anisotropic. In order to be able to compare the experimental numerical diffusivity with the analytical value, the main orientation of the membranes was aligned with the $y$ -axis. The resulting geometries can be found in Figure 5. Two cases have been considered.

• • Case A: In this case, almost perfect layers have been used.
• • Case B: Here we want to estimate the influence of short circuits and more irregular structures.

The geometrical data for both data are provided in Table 1. The corresponding data for the diffusion constants are given in Table 2. Observe that the diffusion in the lipophilic part is anisotropic. This has been used for the numerical experiments. In contrast to that, the homogenized diffusion constant has been determined by using $d_{2,11}$ , only. So we expect a larger error in the experiments with the domain of case B.Figure 4: Contrast amplified reference domain. The black areas indicate membranes

(a)

(b)

Figure 5: Computational domains for case A (a) and case B (b)

	Case A	Case B
$l\_x[\backslash text\{m\}]$	$4.359\; \backslash times\; 10^\{-7\}$	$4.390\; \backslash times\; 10^\{-7\}$
$l\_y[\backslash text\{m\}]$	$0.9125\; \backslash times\; 10^\{-7\}$	$1.568\; \backslash times\; 10^\{-7\}$
$p\_1$	0.8122	0.8139
$p\_2$	0.1878	0.1861

Table 1: Geometric constants

	Value
$d\_1[\backslash text\{m\}^2\backslash text\{s\}^\{-1\}]$	$1.0\; \backslash times\; 10^\{-14\}$
$d\_\{2,11\}[\backslash text\{m\}^2\backslash text\{s\}^\{-1\}]$	$1.0\; \backslash times\; 10^\{-12\}$
$d\_\{2,22\}[\backslash text\{m\}^2\backslash text\{s\}^\{-1\}]$	$1.0\; \backslash times\; 10^\{-10\}$
$K\_p$	$1.26\; \backslash times\; 10^\{-2\}$

Table 2: Diffusion constants

The experimental determination of the effective diffusivity according to Section 3 can be carried out using two approaches:

1. 1. Use the original equation (3) subject to the inner boundary conditions (4) and (5).
2. 2. Use the transformed problem (7) without any inner boundary conditions.

In order to be as close as possible to the original problem we have chosen the first alternative for our experiments. Note that, in the case of $K_p = 1$ , both approaches are identical.

Unfortunately, it is not possible to formulate the inner transmission conditions (4), (5) directly in Comsol Multiphysics. Instead, both conditions have been coupled by a penalty approach as suggested in [3]. For a suitably chosen constant $\kappa$ , (4), (5) is replaced by

$\kappa$ acts as a mass transfer coefficient.

For comparison purposes, even the homogenized problem (10) has been implemented in Comsol Multiphysics.

4.2 Results

The experiments have been carried out using the values

$$\kappa ={10}^{-4},M={10}^{-7},c_0=1,c_1=0.\backslash kappa\; =\; 10^\{-4\},\; \backslash quad\; M\; =\; 10^\{-7\},\; \backslash quad\; c\_0\; =\; 1,\; \backslash quad\; c\_1\; =\; 0.$$

The penalty parameter has been chosen such that both sides of the equations (18) are somehow in balance. The value of $M$ has been chosen such that the outflow has the order magnitude $0.3c_0$ . In case 1, $K_p = 1$ while, in case 2, $K_p = 0.0126$ . The results are summarized in Table 3.

The effective diffusivities given above refer to the steady state. In order to get a feeling for the influence of the homogenization on the transient behavior, we compared the time history of the mean flux out of the domain at the left boundary between the original equation (3) and its

case	hom. constant	exper. constant	rel. difference
1A	1.2284	1.3131	6.9%
1B	1.2258	1.3590	10.9%
2A	1.2312	1.2910	4.9%
2B	1.2286	1.4680	19.5%

Table 3: Homogenized and experimental effective diffusivities scaled by $10^{-14}$ Figure 6: Comparison of the flux at the outflow boundary for the homogenized model (line) and the detailed model (dashed line)

homogenized counterpart (10). For that, the boundary value problem has been solved as before using the initial condition

$$u(x,y)={10}^{-6}\exp (-1000{\left(\frac{x}{2.179\times {10}^{-7}}\right)}^2)\text{at }t=0.u(x,\; y)\; =\; 10^\{-6\}\; \backslash exp\; \backslash left(\; -1000\; \backslash left(\; \backslash frac\{x\}\{2.179\; \backslash times\; 10^\{-7\}\}\; \backslash right)^2\; \backslash right)\; \backslash quad\; \backslash text\{at\; \}\; t\; =\; 0.$$

The value of the experimental effective diffusion has been used in the homogenized problem. The results for the four different cases differ only marginally. As expected from the experiments, the largest differences occur in case 2B. This is shown in Figure 6.

Summarizing, the following sources of errors occur:- • The sizes of the sub-structures are not infinitesimal small;

• • The membrane layers are not ideal planes;
• • In the geometry case B, the membranes are touching the outflow boundary;
• • For computing the homogenized diffusion coefficient, only the normal part of the membrane diffusion tensor has been used;
• • The partition coefficients are handled by a penalty approach.

The size of the numerical errors is negligible compared to the ones given above.

5 Experimental Effective Diffusivities in Random Media

5.1 The Experimental Set-up

The idea for estimating the effective diffusivity in the present case is to use a Monte Carlo simulation. For that, the test domain $D$ of Section 3 is chosen to be the unit cube, $D = (0, 1)^3$ . Let

$$Q=\text{diag}(d_{11},d_{22},d_{33})Q\; =\; \backslash text\{diag\}(d\_\{11\},\; d\_\{22\},\; d\_\{33\})$$

be a fixed diffusion tensor. For a given $N \in \mathbb{N}$ , this cube is subdivided into $N^3$ sub-cubes

$$D_{ijk}=(x_{i-1},x_i)\times (y_{j-1},y_j)\times (z_{k-1},z_k)D\_\{ijk\}\; =\; (x\_\{i-1\},\; x\_i)\; \backslash times\; (y\_\{j-1\},\; y\_j)\; \backslash times\; (z\_\{k-1\},\; z\_k)$$

with $x_i = y_i = z_i = ih$ and $h = 1/N$ . $h$ plays the rôle of $\varepsilon$ in Theorem 2.2. One experiment consists of choosing a realization $d^\varepsilon$ such that

$$d^{\epsilon }{\mathrm{\mid }}_{D_{ijk}}=T_{ijk}Q{T}_{ijk}^{T}d^\backslash varepsilon|\_\{D\_\{ijk\}\}\; =\; T\_\{ijk\}\; Q\; T\_\{ijk\}^T$$

where $T_{ijk} \in SO(3)$ are drawn uniformly distributed in $SO(3)$ .

In order to describe the orientation we will use the Euler angles. Any rotation in $SO(3)$ can be described by three angles, the so-called Euler angles. We will use the convention to first rotate around the $x_3$ -axis by the angle $\alpha$ , then around the (new) $x_1$ -axis by $\beta$ , and finally around the new $x_3$ -axis by $\gamma$ . This can be described formally by

$$T=R_3(\gamma )R_1(\beta )R_3(\alpha ),\alpha ,\gamma \in (0,2\pi ),\beta \in (0,\pi ),(19)T\; =\; R\_3(\backslash gamma)R\_1(\backslash beta)R\_3(\backslash alpha),\; \backslash quad\; \backslash alpha,\; \backslash gamma\; \backslash in\; (0,\; 2\backslash pi),\; \backslash quad\; \backslash beta\; \backslash in\; (0,\; \backslash pi),\; \backslash quad\; (19)$$

where

Let $\mu$ denote the Haar measure on $SO(3)$ . Its density has the simple form

$$d\mu =\frac{1}{8{\pi }^2}\sin \beta d\alpha d\beta d\gamma d\backslash mu\; =\; \backslash frac\{1\}\{8\backslash pi^2\}\; \backslash sin\; \backslash beta\; d\backslash alpha\; d\backslash beta\; d\backslash gamma$$with respect to the Lebesgue measure on $(0, 2\pi) \times (0, \pi) \times (0, 2\pi)$ .

This way, the expectation value of $d_{\text{eff}}^\varepsilon$ can be estimated for given $N(\varepsilon)$ . The computation of $u_{\text{out}}$ and $N_{\text{average}}$ consists essentially of the evaluation of integrals

$${\int }_{{\mathrm{\Gamma }}_1}ud\mathrm{\Gamma }\backslash int\_\{\backslash Gamma\_1\}\; u\; d\backslash Gamma$$

for functions $u \in H^1(G)$ . Since this is a continuous linear functional, we obtain

$$d_{\text{eff}}^{\epsilon }⟶d_{\text{eff}}\text{ for }\epsilon ⟶0d\_\{\backslash text\{eff\}\}^\backslash varepsilon\; \backslash longrightarrow\; d\_\{\backslash text\{eff\}\}\; \backslash text\{\; for\; \}\; \backslash varepsilon\; \backslash longrightarrow\; 0$$

by using (17). Since there are no preferred directions in this setting, the diffusion is isotropic such that $d_{\text{eff}}$ is a scalar.

5.2 Results in 2D

In the two-dimensional setting, an analytic solution of the random homogenization problem is known. Let

$$Q=\text{diag}(d_{11},d_{22})\mathrm{.}Q\; =\; \backslash text\{diag\}(d\_\{11\},\; d\_\{22\}).$$

The effective diffusivity is the scalar [5, p. 235]

$$d_{\text{eff}}=(d_{11}{d}_{22}{)}^{1\mathrm{/}2}\mathrm{.}d\_\{\backslash text\{eff\}\}\; =\; (d\_\{11\}d\_\{22\})^\{1/2\}.$$

We will carry out the experiment described above in the two-dimensional setting in order to obtain a certain gauge for its three-dimensional equivalent.

The two-dimensional counterpart of the experiment described in Section 5.1 is to choose $D = (0, 1)^2$ which will be subdivided, for a given $N \in \mathbb{N}$ , into sub-squares

$$D_{ij}=(x_{i-1},x_i)\times (y_{j-1},y_j)D\_\{ij\}\; =\; (x\_\{i-1\},\; x\_i)\; \backslash times\; (y\_\{j-1\},\; y\_j)$$

with $x_i = y_i = ih$ and $h = 1/N$ . The realizations $d^\varepsilon$ are now described by

$$d^{\epsilon }{\mathrm{\mid }}_{D_{ij}}=T_{ij}Q{T}_{ij}^{T}d^\backslash varepsilon|\_\{D\_\{ij\}\}\; =\; T\_\{ij\}\; Q\; T\_\{ij\}^T$$

where $T_{ij} \in SO(2)$ are sampled uniformly distributed in $SO(2)$ . The elements of $SO(2)$ are simple rotations described uniquely by an angle $\varphi \in [0, 2\pi)$ ,

The Haar measure $\mu$ on $SO(2)$ has the density $d\mu = \frac{1}{2\pi} d\varphi$ with respect to the Lebesgue measure on $(0, 2\pi)$ .

The experimental results for

$$d_{11}=1,d_{22}=10,d_{\text{eff}}=3.1623d\_\{11\}\; =\; 1,\; \backslash quad\; d\_\{22\}\; =\; 10,\; \backslash quad\; d\_\{\backslash text\{eff\}\}\; =\; 3.1623$$

are provided in Table 4. We can draw the following conclusions:- • The main parameter for the accuracy of the estimation of the effective diffusivity is $N$ . This isn't hardly surprising.

• • For a given $N$ , the sample size has only a minor influence on the accuracy. Once a certain number of trials has been reached the accuracy does not become better. The optimal sample size seems to be independent of $N$ .
• • The standard deviation for sufficiently large sample sizes roughly halves while doubling $N$ . This indicates a linear rate of convergence.
• • In all experiments, the mean value of the experimental effective diffusivity is an overestimation of the exact value.
• • If the sample size is too small, the standard deviation is misleading small.
• • The experiments in Section 4 indicated an error in the order of magnitude of 5% – 10% between the theoretical homogenized diffusivity and the experimentally observed. These results suggests to use a value of $N = 20$ and a sample size of at least 15 trials.

5.3 Results in 3D

Finally, the experimental estimation in the three-dimensional case has been carried out. Unfortunately, the geometry handling in Comsol Multiphysics has led to a severe restriction on how large $N$ can be. Although the machine used had enough memory installed (16 GB), the Java heap space got exhausted rather soon. Moreover, the geometry analysis was surprisingly time-consuming compared to the assembly and solution process.

For the experiments, we have chosen the diffusion coefficients

$$d_{11}=9,d_{22}=25,d_{33}=1.d\_\{11\}\; =\; 9,\; \backslash quad\; d\_\{22\}\; =\; 25,\; \backslash quad\; d\_\{33\}\; =\; 1.$$

Note that an analytical solution of the homogenization problem is not known. The results are given in Table 5.

6 Conclusions

The present paper explains the homogenization strategy which has been used to derive effective equations for modelling the detailed metabolism in mammalian cells. The cytoplasm has been modelled assuming that three different length scales can be observed. For going from the smallest to the medium scale, an analytic homogenization technique is used. By comparing the analytic effective diffusion constant with results from numerical simulations on real cell geometries taken from photographs an error of 5% – 20% has been observed. Given the accuracy of the known diffusion constants in the lipophilic and aqueous parts of the cytoplasm this accuracy appears to be sufficient.

$N$	sample size	mean	standard deviation	abs. error
20	5	3.1693	0.0788	0.0070
	10	3.2689	0.1604	0.1066
	15	3.1834	0.1448	0.0211
	30	3.2225	0.1294	0.0602
	60	3.2059	0.1569	0.0436
	90	3.1973	0.1448	0.0350
	120	3.1708	0.1516	0.0085
	150	3.2109	0.1371	0.0486
	180	3.1946	0.1431	0.0323
	200	3.1971	0.1451	0.0348
40	5	3.2377	0.0672	0.0754
	10	3.2272	0.0602	0.0649
	15	3.2343	0.0675	0.0720
	30	3.2380	0.0789	0.0757
	60	3.2496	0.0722	0.0873
	90	3.1907	0.0746	0.0284
60	5	3.1950	0.0276	0.0327
	10	3.1870	0.0487	0.0247
	15	3.1896	0.0420	0.0273
	30	3.1916	0.0417	0.0293
80	15	3.2009	0.0305	0.0386

Table 4: Experimental effective diffusivities in 2D for $d_{11} = 1$ , $d_{22} = 10$ , $d_{\text{eff}} = 3.1623$

$N$	sample size	mean	standard deviation
4	5	7.6753	0.4767
	10	7.2391	0.6292
	15	7.3391	0.6667
	20	7.5785	0.8431
	30	7.5144	0.7630
8	5	8.1298	0.1226
	10	8.1251	0.2088
	15	8.0147	0.2914
	20	8.0910	0.2395
	30	8.0490	0.2193
10	5	8.3499	0.1448
	10	8.2605	0.1769
	15	8.2741	0.1523
	20	8.2783	0.1522
	30	8.3131	0.1524
16	5	8.6457	0.0943
	10	8.7546	0.1005
	15	8.6834	0.0748
	20	8.6453	0.0845
	30	8.6787	0.0752
20	5	8.7419	0.1162
	10	8.7383	0.0622
	15	8.7214	0.0616
	20	8.7412	0.0505
	30	8.7281	0.0596

Table 5: Experimental effective diffusivities in 3D for $d_{11} = 9$ , $d_{22} = 25$ , $d_{33} = 1$For the step from the medium scale to the large scale, a random homogenization technique has been used. Mathematically the effective diffusivity is known to exist. In the present paper an algorithm has been developed and tested for estimating the homogenized diffusion constant on the large scale. However, the computation times in Comsol Multiphysics have become very large (up to one week on a compute server based on a 2GHz AMD Opteron processor) for a reasonable setup such that alternative solution techniques should be investigated.

The critical assumption in the last step is that about the probability distribution of the structures on the intermediate scale. Its validity can probably only be justified by comparison to biochemical experiments.

More detailed results can be found in [1].

Acknowledgement. The authors are grateful to Dr. Holger Jastrow, Mainz, for providing us with a high-resolution photograph of cell organelles.

References

• [1] Marry Chriz Cabauatan-Villanueva. On the computation of effective diffusivities in the cytoplasm of a cell. Master's thesis, Chalmers University of Technology, Göteborg, Sweden, 2007.
• [2] COMSOL AB, Stockholm. Comsol Multiphysics 3.3, 2007.
• [3] COMSOL AB, Stockholm. Comsol Multiphysics 3.3, Chemical Engineering Module, 2007.
• [4] Kristian Dreij, Ralf Morgenstern, Bengt Jernström, and Michael Hanke. A Homogenization Method for Efficient Calculation of Diffusion and Reactions of Lipophilic Compounds in Complex Cell Geometry. In preparation.
• [5] V.V. Jikov, S.M. Kozlov, and O.A. Oleinik. Homogenization of differential operators and integral functionals. Springer-Verlag, Berlin, 1994.
• [6] Lars-Erik Persson, Leif Persson, Nils Svanstedt, and John Wyller. The homogenization method. Studentlitteratur, Lund, 1993.
• [7] Leif Persson. Computing effective thermal conductivities of composite materials by the homogenization method. Examensarbete, Tekniska Högskolan, Lund, 1986.