MSc Physics and Astronomy
Track Theoretical Physics
Master Thesis
Numerical Analysis of the Influx of Diffusive Moving Particles
into Partially Absorbing Non-Spherical Cells
by
Felicio Gordilho Fernandes
5822718
May 2021
60 ECTS
Research conducted between the 1st of September 2015 and the 28th of February 2017
Supervisor/Examiner:
Examiner:
Prof. Dr. Bela Mulder
Dr. Greg Stephens
AMOLF
Living Matter Program
Numerical Analysis of the Influx of Diffusive
Moving Particles into Partially Absorbing
Felicio Gordilho Fernandes ALL RIGHTS RESERVED
Figure 1.1on page 5 from Cell Biology by the Numbers [11], illustrated by Nigel Orme.
Reprinted with permission.
Felicio Gordilho Fernandes
Numerical Analysis of the Influx of
Diffusive Moving Particles into Partially
Abstract
Biological cells can use diffusion to take up material at almost no energy cost. In some of these processes, material can only enter at specific transfer sites. In this thesis, I investigate the effects of the shape of rod-like shaped cells and the distribu-tion of transfer sites on such uptake processes.
I developed a computational method to sample the stochastic absorption process yielding the steady-state rate of the influx of the particles. The algorithm combines three techniques for different regions in space to simulate the behaviour of the dif-fusing particles. For the absorption of the particles, I represent the distribution of transfer sites by a local uptake rate, which in turn is linked to the local absorp-tion probability. My results match the analytical predicabsorp-tions for the spherical cells, though with a small systematic error.
I tested different absorption probabilities proposed in the literature. I find that the absorption probability by Singer et al. (2008) has the best fitting results. For the rod-like shaped cells, I find that their morphology does play a role in the influx of particles. The shape, however, is considerably more important in this than the dis-tribution of local uptake rates. The exact relation between these two characteristics and the influx of particles is non-trivial and remains to be further explored.
Acknowledgements
I am grateful to AMOLF and the Dutch Research Council (NWO) for making my internship possible.
I thank my supervisor Bela Mulder for the opportunity to join his warm group. Your mentorship, involvement, and feedback have helped me throughout this project and beyond.
I thank the rest of the Mulder group, all of the ten Wolde group, and the people of AMOLF in general, for all the discussions, work-related and non-work related. You brought a great atmosphere to my stay at AMOLF. I specifically thank Marco Saltini and Harmen Wierenga for all their input.
I thank all my family, friends, and peers who have supported me in the completion of writing this thesis and in life.
Contents
1. Introduction 1
1.1. Diffusion is a distributive mechanism . . . 1
1.1.1. Brownian motion . . . 1
1.1.2. Diffusion equation . . . 2
1.2. Cells interact with diffusion processes . . . 3
1.2.1. Membrane. . . 3
1.2.2. Interactions . . . 4
1.2.3. Possible adaptive trait: shape . . . 5
1.2.4. Possible adaptive trait: density of transfer sites . . . 6
2. Methodology 7 2.1. Mathematical description . . . 7
2.1.1. Representing free diffusion. . . 7
2.1.2. Boundary conditions characterise other systems . . . 8
Absorbing boundary condition (Dirichlet boundary condition). 8 Reflecting boundary condition (Neumann boundary condition) 9 2.1.3. Formulaic expression of the problem (Robin boundary condition) 10 2.1.4. Analytical method insufficient. . . 11
2.2. Computational methods . . . 12
2.2.1. Finite element method . . . 12
2.2.2. Monte Carlo method . . . 12
A lattice random walk or a Gaussian random walk . . . 12
2.2.3. Simulation scheme of this research . . . 13
Surrounding sphere. . . 14
Inside Gaussian layer. . . 15
Outside Gaussian layer. . . 16
Outside surrounding sphere . . . 17
Complete layout . . . 18
2.3. Statistical methods . . . 19
3. Results 21 3.1. Validation algorithm . . . 21
3.1.1. Thickness Gaussian layer . . . 26
3.1.2. Comparison with alternative proposed absorption probabilities 28 3.2. Simulations spherocylinder . . . 32
3.3. Discussion . . . 33
3.3.1. Top performance absorption probability equation . . . 33
3.3.2. Discrepancy results and analytic solution . . . 34
3.3.3. Low estimates standard error of the mean . . . 34
3.4. Outlook and future research . . . 35
4. Conclusion 37
B.3. Absorbing boundary conditions . . . 42
B.3.1. Steady-state solutions . . . 43
B.3.2. Time-dependent solutions . . . 43
B.3.3. Implementation . . . 47
Cut-off infinite sum . . . 48
B.4. Reflecting boundary conditions . . . 48
B.4.1. Steady-state solutions . . . 48
B.4.2. Time-dependent solutions . . . 49
B.4.3. Implementation . . . 52
Cutoff infinite sum . . . 52
Appendix C. One-dimensional Monte Carlo simulations 53 C.1. Quasi-continuous: Gaussian jumps . . . 53
C.2. Pure discrete: Random walker. . . 54
C.3. The simulated systems and their analytical solutions . . . 55
C.4. The simulations. . . 56
Appendix D. Stochastic properties of a one-dimensional random walker 59 Appendix E. Units of the local absorption rate 63 Appendix F. Steady-state solution sphere Robin B.C. 65 F.1. Steady-state solution sphere Robin B.C. . . 65
F.2. Surface flux . . . 66
Appendix G. Spherocylinder 69 G.1. Distribution absorption rate . . . 69
G.1.1. Assumption one . . . 70
G.1.2. Assumption two . . . 72
G.1.3. Implementation . . . 73
Appendix H. Simulating the prolate ellipsoid 75 H.1. Fukushima method of prolate ellipsoid (a = b < c = 1) in normal vector coordinates . . . 75
H.2. Fukushima method normal vector coordinates . . . 79
H.3. Approx of θ near zi = 0 for prolate ellipsoid (a = b < c = 1) . . . 86
H.4. Exact analytical solution for tan(θi) when ρi = 0 for prolate ellipsoid (a = b < c = 1) . . . 88
1. Introduction
This thesis aims at investigating whether cells that depend on specific diffusion processes can be optimized for these processes through their morphology. First I discuss diffusion itself as a concept and how it can be described by formulas. Then I will look into how cells protect themselves against diffusion and how they can use it, so that I can formulate a hypothesis.
1.1 Diffusion is a distributive mechanism
Diffusion is the emergent phenomenon that a concentration of a substance in a medium will disperse equally over the medium, given enough time and no physical restraints, as a result of the non-directional movement of the separate indistinguish-able particles. This can be illustrated by a simple example. Imagine two gasses, so free moving molecules, already at thermal equilibrium in two adjacent compart-ments only separated by a slide door. Once the slide door is opened, there will not be advection — due to the thermal equilibrium — yet the two gasses will diffuse into one another and become a homogeneous mixture over time by cause of the random motion of their molecules.
Diffusion will not always occur when different substances come into contact with each other. It is necessary that they can mix — without adding energy into the sys-tem externally. An example for substances that do not mix is when one substance is hydrophobic and the other one is hydrophilic. This is due their polarity.
1.1.1 Brownian motion
Brownian motion is the underlying mechanism that causes diffusion. Particles in-trinsically have at least some kinetic energy. This means that particles are funda-mentally not completely stationary. In classical mechanics, this kinetic energy is as-sociated to temperature through thermal fluctuations. Thermal fluctuations are the random changes in possible states of a system in equilibrium. As the temperature increases, the thermal fluctuations become larger.
When particles are in a medium in which they can move, they can hit other particles and change their direction. This can make the movement of particles erratic. This phenomenon is called Brownian motion, after the Scottish botanist Robert Brown (1773–1858) [1], and was first explained by Albert Einstein in 1905 [2,3].
Brownian motion causes particles — when they are not strongly bound to some-thing — to disperse over a medium, forced by all the collisions they endure.
with j the diffusive flux of particles: the amount of diffusing material passing through a unit area in a unit of time; and c the concentration of particles: the number of particles per volume2. The continuity equation (1.1) expresses that the change in
time in the concentration of particles, c, is equal to the change in space in the flux of particles, j. Now one can relate this flux of particles j to diffusion via Fick’s first law:
j =−D∇c, (1.2)
which relates the flux of particles, j, to the concentration gradient,∇c, via a pro-portionality factor: the diffusion coefficient D.
Therefore, the diffusion equation looks in its most general form like: ∂c(r, t) ∂t =∇ · D c(r, t), r∇c(r, t) , (1.3)
with c(r, t) the local — at r — concentration of particles of the substance at time t, and D c(r, t), r the diffusion coefficient. If D does not depend on the local con-centration of particles, equation (1.3) is a linear equation. If D is constant — then also called the diffusion constant — the medium is homogeneous and the diffusion equation (1.3) becomes:
∂c(r, t)
∂t = D∇
2c(r, t). (1.4)
In three-dimensional Cartesian coordinates the vector differential operator∇ is:
∇ = ∂ ∂r = ˆ ex ∂ ∂x+ ˆey ∂ ∂y + ˆez ∂ ∂z , (1.5)
and hence the diffusion equation (1.4) becomes: ∂c(r, t)
∂t = D
∂2
∂r2c(r, t). (1.6)
Since the three Cartesian dimensions are orthogonal and c(r, t) is uniform in space for free diffusion, a solution for one-dimensional diffusion can easily be expanded to three dimensions by multiplying the dimensional part of the one-dimensional solution for all independent dimensions (see section A.2on page 40).
It is possible to nondimensionalize the one-dimensional diffusion equation (A.2) for finite systems (see page 41, equations (B.1) and (B.2)); the characteristic units to do this are the diffusion constant and the length of the system. By nondimen-sionalizing, the diffusion constant can literally be taken out of the equation, and
1In case of a source (material generated): σ > 0, and in case of a sink (material removed): σ < 0.
2Molar density for the concentration would work too: the number of moles per volume; in the
1.2. Cells interact with diffusion processes 3
one can easily compare different systems with the same mechanics. The general dynamics in a nondimensionalized system are described by the appropriate diffusion equation (B.3):
∂c(ξ, τ )
∂τ =
∂2c(ξ, τ )
∂ξ2 , (1.7)
with ξ the dimensionless position parameter and τ the dimensionless time para-meter, and boundaries at ξ = 0 and ξ = 1. In appendixB,Solution 1D nondi-mensionalized diffusion, on page41 this is used to find solutions for a system with absorbing boundaries, and one with reflecting boundaries.
1.2 Cells interact with diffusion processes
Biological cells exist at a scale at which molecular diffusion occurs. To give an idea, for a diffusive moving protein of about a nanometre in length, in a cell with a water-like interior, so a viscosity η of about 9.58· 10−4 Pa· s at an absolute temperature T of 295 K, and approximating the protein as a spherical particle, we have with the Stokes–Einstein equation for the diffusion constant:
D = kBT 6πηr = 1.380649· 10 −23· 295 6π· 9.58 · 10−4· 1 · 10−9 m 2· s−1 = 2.25548 . . .· 10−10m2· s−1 = 2· 10−10 m2· s−1, (1.8)
with kB the Boltzmann constant, and r the radius of the spherical approximation of the protein. The time scale to get from two points in the cell a micrometre apart would be (for three-dimensional diffusion):
∆t = R 2 6D = 1· 10 −92 6· 2.25548 . . . · 10−10 s = 7.38942 . . .· 10−10 s = .7 ms, (1.9)
with R the distance between the two points.
Because biological cells exist at such a scale, they therefore interact with the phe-nomena related to diffusion. As a consequence, phephe-nomena related to biological cells can be studied with diffusion in mind. That is to say, we can question to what extent a particular biological phenomenon is (mainly) influenced by diffusion.
1.2.1 Membrane
The membrane is the outer layer of the cell. It generally prevents certain particles from diffusing into or out of the cell. Without it, the contents of a cell could simply diffuse away.
The membrane is a lipid bilayer. This means that it is a layer of two amphipathic particles opposing each other. So these are two particles with each a hydrophilic
protects the cell in another way by keeping its contents from getting into contact with just any other particles in the medium.
1.2.2 Interactions
Cells do need to move particles in or out however, as they need to take up nutrients and get rid of waste. Diffusion plays a role in the movement of the particles around a cell and can therefore be of importance in these interactions.
There are many ways different cells can transfer material across the bilayer. I will focus on the simplest one, being passive transport. In passive transport, the cell does not use its energy to drive the transport. Instead, the transport is driven by diffusion: the net flux of particles from a high concentration to a low concentration so that these concentrations become more alike.3 For different particles, the high
concentration can be either in or outside the cell.
When particles are small enough and do not have (much) charge, they can diffuse through the cell membrane by themselves. This process is called simple diffusion. Larger or charged particles need help to cross the bilayer for passive transport to happen.4 This is done by specific molecules in the membrane that function as trans-fer sites. This process is known as facilitated diffusion. There are generally two types of transfer sites through which this can take place.
Firstly, it can take place via channel proteins in the bilayer. Channel proteins allow specific particles to tunnel through the membrane, for example porins in the mem-branes of bacteria and eukaryotes [4, p. 635]. Porins are filled with water and allow the transport of some small hydrophilic particles. Narrow porins can be accessible to certain molecules only, such as aquaporins that facilitate osmosis [4, pp. 673–675]. Osmosis is the transport of water across the membrane from high to low concentra-tion of water, so from low to high concentraconcentra-tion of solute. Aquaporins are so narrow that water molecules have to pass in a single line, and hydrated ions cannot pass at all. For the ions, there are ion channels [4, pp. 667–693]. Which are selective to spe-cific ions and they do not have to be permanently open. This way they can regulate a concentration. Examples of these ions are sodium: Na+, potassium: K+, calcium: Ca2+, and chloride: Cl−. Because of the charge of these ions and the ability to regulate the concentration, ion channels play a role in the electrical transmission necessary for nerve and muscle cells.
Secondly, facilitated diffusion can also take place with the help of carrier proteins in the bilayer. Carrier proteins bind to specific particles, change shape, and move the particles to the other side of the membrane. Examples of the use of carrier proteins are some transport mechanisms of sugars across cell membranes [5,6]. Oxygen
3This is different from active transport, which uses energy to fuel transport from a lower
concentration to a higher concentration.
4Some particles can diffuse through the membrane by themselves in principle, but their transfer
1.2. Cells interact with diffusion processes 5
Figure 1.1.: Examples of different cell shapes [11, p. 51], reprinted with permission. The original caption reads: “A gallery of microbial cell shapes. These drawings are based upon microscopy images from the original literature. (A) Stella strain IFAM1312 (380); (B) Microcyclus (a genus since
renamed Ancylobacter) flavus (367); (C) Bifidobacterium bifidum; (D) Clostridium cocleatum; (E) Aquaspirillum autotrophicum; (F) Pyroditium abyssi (380); (G) Escherichia coli; (H) Bifidobacterium sp.; (I) transverse section of ratoon stunt-associated bacterium; (J) Planctomyces sp. (133); (K) Nocardia opaca; (L) Chain of ratoon stunt-associated bacteria; (M) Caulobacter sp. (380); (N) Spirochaeta halophila; (O) Prosthecobacter fusiformis; (P) Methanogenium
cariaci; (Q) Arthrobacter globiformis growth cycle; (R) gram-negative Alphaproteobacteria from marine sponges (240); (S) Ancalomicrobium sp. (380); (T) Nevskia ramosa (133); (U) Rhodomicrobium vanniellii; (V) Streptomyces sp.; (W) Caryophanon latum; (X) Calothrix sp. (Y) A
schematic of part of the giant bacterium Thiomargarita namibiensis (290). All images are drawn to the same scale. (Adapted from K. D. Young, Microbiology & Molecular Bio. Rev., 70:660, 2006.)” [11, p. 51].
and carbon monoxide also use carriers: haemoglobin in blood and myoglobin in muscles [7–10].
1.2.3 Possible adaptive trait: shape
A possible way for cells to be optimized for diffusive processes could be the shape that they have. Cells can be found in a large variety of shapes, for example spher-ical (such as coccus bacteria), rod-like (bacilli), or more ellipsoid cells (coccobacilli). More examples of cell shapes are shown in figure1.1 [11, p. 51] on page 5.
It is known that many prokaryotic and eukaryotic cells have optimal sizes for in-tracellular diffusion with respect to the number of macromolecules that need to be
∂c(r, t)
∂t =∇ · D(n, r)∇c(r, t)
, (1.10)
one can spot the similarities with Poisson’s equation: − ρ
ε =−∇ · E = ∇ · (∇V ), (1.11)
where ρ is the total volume density, ε the permittivity, E the electric field, and V the electric potential.
Likewise, one can see the resemblance between the equation for the strength of an electrostatic field:
E =−∇V, (1.12)
and Fick’s first law, that was used on page2:
j =−D∇c. (1.13)
Following these comparisons, one can construct an analogy between a conductor in an electrostatic field and a cell in a diffusive medium: at places where the curvature of the cell is higher, the influx of material will be higher too.5
1.2.4 Possible adaptive trait: density of transfer sites
Cells could hypothetically be optimized for facilitated diffusion by having a high density of uptake channels or carrier proteins at certain locations. One can assume that when a particle requires facilitated diffusion to get into a cell, the probability of accessing the cell will increase if there are more transfer sites close to the nearing particle. Hence, I hypothesize that if the curvature of a cell indeed benefits the influx of particles, having more transfer sites at the more curved areas should be beneficial compared to a uniform distribution.
5It should be noted that the strength of the electrostatic field is caused by the charge density
on the conductor. So this is where the analogy with the cell in a diffusive medium breaks. As an aside, the analytic relation between the charge density at a point of a conductor and the curvature, is still an object of study, see for example: the study On the dependence of charge density on
surface curvature of an isolated conductor from 2016 [13]. To add to that, what the ideal shape is for a lightning rod has been under debate for a long time. See for example Nikola Tesla’s patent
Lightning-protector from 1918 [14], and the more recent papers The case for using blunt-tipped
lightning rods as strike receptors [15] and On the optimum rod geometry for practical lightning
2. Methodology
The research questions posed in the previous chapter could be answered if one could measure how much material a cell could take up. Cells with different traits could then be compared. One could measure how much material a rod-like shaped cell takes up when it has more transfer sites on its tips than on the centre of its body. Then, it could be compared that with the taking in of a complete symmetric spher-ical cell with a uniform distribution of transfer sites.
Conducting these measurements would be time-consuming and require a complic-ated experimental set-up, including the right microscopes and specimens. Fortu-nately, there is a different, more time and cost-economic option.
2.1 Mathematical description
I can articulate my biological system in mathematical terms. In this research, I simplify the intricate process of facilitated diffusion by disregarding the complex dynamics associated with the processes at the transfer sites. Instead, I state that when a particle hits the membrane, there is a certain probability, Pabs, that this particle will be absorbed. If the particle is not absorbed, it will be reflected with probability Pref = 1− Pabs. I need to be able to let this Pabs change over the surface of the cell, correlating to the local number of transfer sites. I will represent this sought behaviour by a boundary condition on my system so that I can work with it in a physics’ framework. In order to do this I need an understanding of different simple models first: of a freely diffusive system, of one with absorbing boundaries, and of one with reflecting boundaries.
2.1.1 Representing free diffusion
With a solution for free diffusion I can describe the general dynamics of a diffusive moving particle. This solution for free diffusion with a constant diffusion coefficient is discussed in appendixA (page 39). From this procedure for c(r, t) I know that the concentration profile of a collection of particles in a one-dimensional1, x = (x), free diffusive system, with at t = 0 all particles at x = 0, is described by a Gaussian function: c(x, t) = N√ 1 2πσ2e −x2 2σ2, (2.1) with σ2 = 2Dt, (2.2)
and N the number of those particles (see equations (A.3) and (A.4)).
1The one-dimensional position is written here explicitly as a vector to make its relation to the
three-dimensional case more natural. If I am truly consistent and pedantic I express time in its one-dimensional vector space too: t = (t), through which the concentration profile is denoted as c(x, t). But this would be unnecessarily unorthodox.
with σ2 the variance.
As noted on page2, the three-dimensional, r = (x, y, z), generalization is the product of the solutions of the three independent dimensions (see also equations (A.6) and (A.7)): Pparticle(r, t) = 1 (2πσ2)3/2e −x2+y2+z2 2σ2 , (2.4) with σ2 = 2Dt. (2.5)
With these equations I can represent a freely diffusive moving particle, such that with a collection of those, I get diffusion.
2.1.2 Boundary conditions characterise other systems
Additional dynamics to a diffusive system are provided by specific conditions on the boundaries of that system. The two simplest examples are the absorbing and reflect-ing boundary conditions, which are forms of the Dirichlet and Neumann boundary conditions, respectively. Their functioning can be illustrated by applying them to a nondimensionalized one-dimensional system (equation (1.7)). This gives insight into the boundary conditions I need for my problem. Both cases will be mentioned here in short, for all the details see sectionsB.3 and B.4(pages42 and 48).
Absorbing boundary condition (Dirichlet boundary condition)
When the concentration of particles on the boundary is fixed to zero, one has the absorbing boundary condition. This is because when a particle hits the boundary, it is forced to disappear. The absorbing boundary condition is a form of the Dirichlet boundary condition, where the value of a function is fixed on the boundary of the domain:
a· f(r) = g(r) on the boundary of the system. (2.6)
This is the absorbing boundary condition when applied to the concentration of particles and one sets g(r) = 0, because then no particles can exist on the boundary due to the inevitable outflux.
Thus, when I apply this to my nondimensionalized one-dimensional system, I have:
c(0, τ ) = 0 = c(1, τ ). (2.7)
I find that the number of particles, c(ξ, τ ), at a certain position ξ at a certain time τ is (equation (B.47)): c(ξ, τ ) = ∞ X m=1 2 sin(mπξstart)N e−m2π2τsin mπξ, (2.8)
2.1. Mathematical description 9
with ξstart the starting position of all the particles in the system, and N the total number of particles in the system.
The absorbing boundary condition in three dimensions becomes:
c(r, t) = 0, r∈ ∂Ω, (2.9)
with c(r, t) the number of particles and ∂Ω the boundary of some domain. Reflecting boundary condition (Neumann boundary condition)
When the value of the first derivative of the concentration of particles with respect to the position is fixed to zero on the boundary, one has the reflecting boundary condition. This is because the number of particles at the boundary cannot change, so they cannot leave the system at that point. The reflecting boundary condition is a form of the Neumann boundary condition, where the value of the first derivative of a function is fixed on the boundary of the domain:
b·∂f (r)
∂n = g(r) on the boundary of the system, (2.10)
where ∂
∂n is the normal derivative — a derivative taken orthogonal to a surface — at r on the boundary. Equation (2.10) is the reflecting boundary condition when applied to the concentration of particles and one sets g(r) = 0 again, because now the outflux on the boundary is zero and no particles can leave the system.
Hence, when I apply this to my nondimensionalized one-dimensional system, I have: ∂c(ξ, τ ) ∂ξ ξ=0 = 0 = ∂c(ξ, τ ) ∂ξ ξ=1 . (2.11)
I find that the number of particles, c(ξ, τ ), at a certain position ξ at a certain time τ is (equation (B.79)): c(ξ, τ ) = N + ∞ X m=1 2 cos(mπξstart)N e−m2π2τcos mπξ, (2.12)
with, again, ξstart the starting position of all the particles in the system, and N the total number of particles in the system.
For the reflective boundary condition in higher dimensional spaces, the normal de-rivative is of importance. For a function f and normal direction n the normal deriv-ative is:
∂f
∂n =∇f(r) · n = ∇f(r) · ˆn(r), (2.13)
with ˆn(r) the unit normal explicit in point r.
The reflective boundary condition in three dimensions is then:
∇c(r, t) · ˆn(r) = 0, r ∈ ∂Ω. (2.14)
with once more c(r, t) the number of particles and ∂Ω the boundary of some do-main.
I can represent this with the Robin boundary condition : a· f(r) + b · ∂f (r)
∂n = g(r) on the boundary of the system, (2.15) where a and b can, in principle, be functions too. Clearly, the Robin boundary con-dition is a combination of the Dirichlet and Neumann boundary concon-ditions (equa-tions (2.6) and (2.10)), they are extreme cases of the Robin boundary condition. If I set g(r) = 0 on the boundary, I can rewrite the Robin boundary condition to the following form:
∂f (r)
∂n = k(r)f (r), r∈ S, (2.16)
with S the boundary of the system, and k(r) the local Robin parameter on that boundary.
Note that, when k(r) is equal to zero, equation (2.16) is just the reflecting bound-ary condition (2.14); and when 1/k(r) → 0, it becomes the absorbing boundary condition (2.9). Hence, this form of the Robin boundary condition can be seen as a combination of these boundary conditions, with k(r) giving the proportion between them. That way, the Robin boundary condition can be used to describe both re-flections and absorptions of particles. The local Robin parameter will become an indicator of how high the density of transfer sites is. So that when there are many, k(r) is high, and the boundary condition is more absorbing. When there are few transfer sites, k(r) is low, and the boundary condition is more reflecting.
Applying the Robin boundary condition to my biological system, I get that the diffusive flux (equation (1.2)) through my absorbing cell is proportional to k(r):
D∇c(r, t) · ˆn(r) = k(r)c(r, t), r ∈ ∂Ω, (2.17)
with D still the diffusion coefficient, c(r, t) the number of particles, ˆn(r) the unit normal, and ∂Ω the boundary of my cell. The Robin parameter k(r) is now the local absorption rate, with different units than in equation (2.15), namely [k(r)] = m s−1 (see appendix Eon page 63). When the density of transfer sites is high, the local absorption rate will be high, and vice versa.
Now that I have the Robin boundary condition to describe a partially absorbing surface, I can try to find the total number of particles that enter. Analogous to in-tegrating the mass flux over a surface and time interval to find the total amount of mass going through that surface in that period, one can integrate the diffusive flux to find the total number of particles, n, that went through some surface S between
2Named after French mathematician Victor Gustave Robin (1855–1897).
It is also known as the impedance boundary condition, or convective boundary condition. In steady state situations (no time dependence) the Robin boundary condition is similar to the radiation boundary condition. That is exactly the case in my research, so in some of the literature that I use the authors refer to the radiation boundary condition instead.
2.1. Mathematical description 11 t1 and t2: n = ˆ t2 t1 ¨ S j· dA dt =−D ˆ t2 t1 ¨ S ∇c(r, t) · ˆn(r) dA dt. (2.18)
Similarly, the flowing of particles through a surface per unit of time Φj3 is:
Φj,S = ¨ S j· dA =−D ¨ S ∇c(r, t) · ˆn(r) dA. (2.19)
In order to solve this integral, I need a solution for the concentration profile c(r, t). I consider that my cell lives a medium that is at equilibrium, or in steady state, and the concentration profile thus does not change in time:
∂css(r)
∂t = 0. (2.20)
This has the logical effect for the other side of the diffusion equation (1.4):
D∇2css(r) = 0. (2.21)
From here on, the subscript in css(r) is dropped and steady state is implied when the concentration of particles only depends on the position: c(r).
Now I want to solve this equation with the appropriate boundary conditions. Be-cause the cell is small with respect to the system, its uptake is small to the medium at equilibrium, I take that the concentration of particles on the boundary of the do-main that I am considering is unperturbed by this uptake. Hence, I use a Dirichlet boundary condition for the boundary of the total domain:
c(|r| → ∞) = c∞, (2.22)
and I use the Robin boundary condition (2.17) on the cell.
2.1.4 Analytical method insufficient
The problem is that one cannot find a general analytical solution for equation (2.21) with the conditions I use. Only when k(r) in the Robin boundary condition (2.17) is uniform over the surface of a sphere, one can find a solution (see appendixFon page65). That is exactly what I do not desire for my cell. I am interested in a non-uniform proportionality factor, so that it can represent a more absorbing bilayer at certain locations of the membrane. Moreover, I am also curious to shapes different from the sphere. Therefore, I need a different way of solving my problem.
3To denote the integral of a flux j over a surface, sometimes the capital J is used, or I (for
In a finite element analysis, a mesh is created to represent the space with smaller separate pieces, the finite elements. In these elements, the equations are solved. The idea is that because the total mesh approximates the actual studied system, their combined solutions do too.
However, the finite element method has some fundamental (potential) disadvant-ages. In relation to the research problem posed in this thesis, these disadvantages present themselves as the following three drawbacks [17]. The finite element method creates a solution for the complete mesh, thus the complete space. Not all of this information about the complete system is essential however, when you just want to know the total influx into the cell; and therefore this seems to be a misuse of computational power. Moreover, the mesh around the absorbing body is an ap-proximation. A better approximation means a finer mesh, at higher computational costs. Therefore, there is a trade-off between the precision of the results and the minimization of computational power. Next to that, for the outer side of the system, the mesh necessitates a physical boundary that, for my purposes, needs to be at a great distance. This means that the complete space for the simulation has to be large, which is again accompanied by computational costs.
For these reasons, I designed a different method to find out in which cases the cell absorbs the most particles.. Since this method is a Monte Carlo simulation, I will first go into what that is and how one can simulate diffusion; then, I will explain the simulation method of this research in section 2.2.3on page 13.
2.2.2 Monte Carlo method
A Monte Carlo method is a simple simulation scheme in which one simulates the dynamics of only one particle using random sampling, keeps track of the outcome, and repeats. After a sufficient number of simulations, the data can be analysed in the same way as experimentally collected data. This analysis should give insight into the behaviour of the system: the statistical properties of the consequences of the stochastic dynamics inside it.
Just like the results of the finite element method are approximates, so are Monte Carlo simulations in their own way, see section 2.3,Statistical methods, on page 19 on Monte Carlo samples.
As mentioned, to apply the Monte Carlo method, I need to be able to simulate a single particle.
A lattice random walk or a Gaussian random walk
A (one-dimensional) diffusive mover can be simulated with either a discrete model, a random walker on a lattice, or more precisely, with a quasi-continuous model
2.2. Computational methods 13
through making Gaussian jumps.4
In a random walk on a lattice, a particle changes location each time step to a neigh-bouring site at random. The resulting probability to be at a certain location matches the distribution of diffusive moving (after enough steps taken).5 Thus, a collection of random walkers on a lattice can simulate diffusion.
With Gaussian jumps, each time step, a sample is taken from the appropriate Gaus-sian distribution. This sample is added to the current location to determine the new location. This will logically result in the correct distribution, because this cascade of samples will create the correct Gaussian distribution itself.6
These Gaussian jumps, or the Gaussian random walk, will be more precise as it will keep space continuous (up to computer precision). Therefore, it does not need to discretize space as with a random walk on a lattice.
For both types of simulations, however, the resolution of the results depends on the size of the time steps that are simulated. This is because the variance in the position of a particle depends on the time (see equations (A.4), (C.2), and (C.5), and their sections).
The one-dimensional models can be developed into three-dimensional models. As the dimensions are orthogonal to each other, movement in each dimension can be simulated independently. The expansion of the Gaussian jumps to three dimensions is straightforward: one-dimensional samples can be used for each direction, creating the right three-dimensional distribution with a mean squared displacement of 6D∆t. With a random walker on a lattice, the particle needs to be able to make a jump in all three directions. There are different ways to achieve this, but they all share that the distance between the new possible sites and the original position needs to match the average distance travelled per time interval. One way of doing this, is to expand the one-dimensional lattice in one plane with two others for the per-pendicular planes. Making a jump in all three planes results in eight possible sites: the vertices of a cube centred around the original position. If the one-dimensional lattice spacing of √2D∆t is kept between adjacent sites, then the eight possible sites are at a distance√6D∆t as required.
With a three-dimensional model of a diffusive particle, the Monte Carlo method can be applied. However, this is not the full story on the simulation method of this research.
2.2.3 Simulation scheme of this research
Together with my supervisor, Bela M. Mulder, I realized that I do not need all the information of the particles everywhere in the diffusive system to find the answer I am looking for. This property is exploited in my simulations. I created a scheme in which diffusive movement that is not close to the absorbing body, is simulated by sampling an element from a sphere surrounding the particle. When the particle gets sufficiently close to the cell, it can be absorbed; and when it gets far-removed from the cell, it can leave the simulation with a certain probability. Therefore, the possible fates of the particle can be tracked and this is a Monte Carlo simulation.
4The fundamentals of these methods are elaborated on in appendixC,One-dimensional Monte
Carlo simulations, on page53.
5See for details sectionC.2,Pure discrete: Random walker, on page54. 6See also sectionC.1,Quasi-continuous: Gaussian jumps, on page53.
Figure 2.1.: A cell in the middle with a surrounding sphere around it described by equation (2.23).
In the following sections, I will introduce the surrounding sphere, the Gaussian layer, the walk-on-spheres layer, and the absorption and escape probabilities that are essential in my design. After introducing all the elements, I can present a flow chart and depiction of the algorithm in figure2.5.
In section3.3,Discussion, on page 33 I review the disadvantages of my method, also in relation to the finite element method.
Surrounding sphere
You could picture a sphere that surrounds the cell — centred at the origin — at a radius RS from the centre of the cell:
S2(0, RS) ={rS ∈ R
3 :kr
Sk = RS}, (2.23)
see: figure2.1.
In order to reach the cell, a particle needs to pass through the sphere first. This means that it is sufficient to know the probability for a particle to be absorbed when it starts out from the surface of this sphere:
Φj,cell = Pabs(∂Ω|S2(0, RS)) S2(0,R S)· Φjabs,S2(0,RS), (2.24)
where Φj,cell is the number of particles that enter the cell per unit time, Φj
abs,S2(0,RS)
is the number of particles flowing through the surrounding, full absorbing, sphere per unit time, and Pabs(∂Ω|S2(0, RS))
S2(0,R
S) denotes the average probability to
being absorbed at the boundary of the cell, ∂Ω, after getting there, when currently located at a point r ∈ S2(0, RS). It is the average probability because the point is that I assume this is not given by a uniform probability distribution on the bound-ary. I will call this average probability the probabilistic influx Pj, such that:
Pj = Pabs(∂Ω|S2(0, RS)) S2(0,R S), (2.25)
and the conditions of the initial location are implied.
Now, I reach an important realization: to find Pj I do not need to have an exact
solution for the concentration profile in space (2.21). I am only interested if a particle is absorbed, not when. That is to say, it is not essential to know at what position ri a given particle at time ti is, but only whether a particle, released somewhere on the
2.2. Computational methods 15
d∗
Figure 2.2.: A cell with a Gaussian layer of thickness d∗ around it. In the layer a diffusive moving particle is simulated with Gaussian jumps.
sphere: r0 ∈ S2(0, R
S), is going to be absorbed. Thus, I can focus on measuring the influx of particles numerically, whilst using the given conditions (equations (2.17), (2.21), and (2.22)).
As I do not need to know the location of specific particle is at a specific point in time, I do not need to simulate the diffusive process by sampling from Gaussian distributions most of the time. This more meticulous method is only necessary near my absorbing object, as I cannot predict how the particles will interact with that exactly.
Inside Gaussian layer: Gaussian jumps and absorption probability
For the purpose of simulating the interaction between the particle and the cell, I make a layer around the cell, with thickness d∗ (see: figure 2.2). In this layer I take samples for all three dimensions from the one-dimensional probability distribution equation (2.3), to determine the new position. In doing so, I effectively take a sample from the three-dimensional probability distribution, equation (2.4), as described in the last part of section 2.2.2, on page 13.
Once the new position of the particle is inside the cell, I need to decide whether it is absorbed or not. For this, I need an expression for the probability of absorption Pabs. The expression has to depend on the local absorption rate k(r), as that is the indicator of the density of transfer sites, as noted on page10. Since the Robin boundary condition is much less prevalent in (studied) physics than the Dirichlet or Neumann boundary conditions, its implementation is less straightforward. In this research, I used developments by Singer, Schuss, Osipov, and Holcman; and work done by my supervisor Bela M. Mulder.
In their 2008 paper, Singer et al. show the relation between “[t]he radiation para-meter κ(x, t) in the d-dimensional Robin boundary condition and the absorbtion [sic] parameter P (x)” [18, p. 4]7:
κ(x, t) = rP (x)pσn(t), x1 = 0, (2.26)
with the constant diffusion matrix σn(t) = nTσ(t)n, and r = 1/ √
π. They note that “The generalization of the multi-dimensional case to domains with curved
boundar-ies […] is not straightforward […]” [18, p. 17]. Mulder shows that this generalization is possible for my absorbing bodies: when it is “a smooth convex bounding surface, provided the local radii of this surface are large compared to the scale of the free diffusion length in a single time step” [17, p. 9]. He argues to set the constant dif-fusion matrix σn(t) to D/∆t. If I now relabel parameters to the ones I use here, I find that the probability of absorption at the boundary of the absorbing object, for a free diffusive moving particle modelled as taking Gaussian jumps per discrete time
Figure 2.3.: A cell with a surrounding sphere (figure 2.1) and a Gaussian layer (figure 2.2), in between these is the walk-on-spheres layer. On the
surrounding sphere and in the walk-on-spheres layer, diffusive movement is simulated by sampling from spheres described by
equation (2.28). An example of such sampling is depicted with dotted lines.
steps and then able to jump into the object, is: Pabs(k) = k·
r π∆t
D , (2.27)
with k the local absorption rate on the boundary, ∆t the size of the simulated time steps, and D the diffusion coefficient.
If the particle is not absorbed, it will be reflected back into the layer, perpendicular to the closest point of the surface.
In section3.1.1,Thickness Gaussian layer, on page26 I look for an optimal value of d∗ for my simulations. In section 3.1.2,Comparison with alternative proposed absorption probabilities, on page 28 I test equation (2.27) against other propositions for the absorption probability.
Outside Gaussian layer: walk-on-spheres
Further away from the cell, I assume that I can approximate the behaviour of the particles with free diffusion, and because of this, I can use an interesting simulation technique. When I have a particle at time ti at position ri, I can look for the closest
point of the cell to that point: rcell,i, and the distance between those two points: RWi=kri− rcell,ik. Then, I can draw a sphere around ri with radius RWi:
S2(ri, RWi) ={rWi ∈ R
3 :kr
Wi− rik = RWi}, (2.28)
such that a part of the surface of this sphere is inside the Gaussian layer, but none of it is inside the absorbing body.
Then I use that the system is ergodic, which means that I know that at some time tj in the future, however far, the particle will be at a point on the surface of that sphere, so that rj ∈ S2(ri, RWi). Since the exact trajectories of specific particles are
2.2. Computational methods 17
not essential to my research question, I can take any random point on the sphere and just place the particle there in my simulation (see figure2.3). This might seem bad practice intuitively, but all I do is: from all the possible particles that could be represented by a particle at ri, I choose one of those that leave the interior of S2(ri, RWi) at a random rj. Since this could be any of the diffusing particles in the
sphere, I still simulate diffusion. This method is called walk-on-spheres, developed in 1956 by Mervin E. Muller [19]. I use it to quickly simulate the non-interacting part of the behaviour of the particle in my system.
Outside surrounding sphere: escape probability
When the new position of the particle is outside the surrounding sphere, S2(0, RS), I have to decide whether the particle returns to the sphere or leaves the system. For this, I follow a method resembling [20], laid out in [17], where an analogy with electrostatics is made once more. Mulder shows that the probability density of not escaping the system satisfies the Poisson’s equation with a boundary condition such that it is similar to the problem of a negative unit charge placed outside a spherical perfect conductor. He then shows that this has as a result that the probability of returning to any point on the sphere is similar to the induced charge density on this conducting surface. For a particle at ri outside the surrounding sphere, such that
krik > RS, I can define (equation (14) in [17]): P ≡ RS
krik
, (2.29)
and I have for the escape probability:
Pesc = 1− P. (2.30)
If the particle does not depart the simulation, it needs to be placed back somewhere on the surrounding sphere (see: figure2.4). To do so, I define a rotated coordinate system X′Y′Z′ with the Z′-axis aligned with ri, and establish the new position in this system first, after which I transform back to the old coordinate system.
The Z′ coordinate can be found by selecting a value for the cosine of the polar angle and then multiplying with the radius of the surrounding sphere. For this I use equa-tion (18) in [17] (which is related to the induced charge at a specific point on the perfect conducting sphere):
cos(θ) = 2ρ
2P (P2+ 1)− 2ρ(P2+ 1)(P− 1) − (P − 1)2
((2ρ− 1)P + 1)2 , (2.31)
with ρ a random sample from the uniform distribution between 0 and 1. Then:
z′= cos(θ)· RS. (2.32)
Now the new X′ and Y′ coordinates will lie on a circle in a plane centred at z′ with radius:
rx′y′ = p
1− cos(θ)2· R
S. (2.33)
If cos(θ) = −1 ∨ 1, then {x′, y′} = {0, 0}; otherwise, I keep picking a set of two random samplesRx′ and Ry′ from the uniform distribution between−1 and 1 till
RS
Figure 2.4.: A cell with a surrounding sphere (figure 2.1) and a particle outside it. The particle can escape the system or return to somewhere on the surrounding sphere.
from just the disc with radius one centred at the origin, and not from the entire square). Then, I scale both accordingly to determine the coordinates:
{x′, y′} = {Rx′,Ry′} · q rx′y′ Rx′2+Ry′2
. (2.34)
Now{x′, y′, z′} ∈ S2(0, RS).
These new coordinates are then transformed back via a transformation defined by ri, so that the new rj = {x, y, z} ∈ S2(0, RS) and the walk-on-spheres method
recommences. Complete layout
Now that we have all pieces, the entire scheme of my event-driven stochastic al-gorithm looks like figure2.5.
So, a particle can leave the system by either being absorbed by the cell or by es-caping to infinity. By counting the number of times a particle is absorbed under certain conditions, I get a data point for the probability to be absorbed under those conditions. This probability is the probabilistic influx Pj mentioned previously,
equation (2.25), and the simulation gives a Monte Carlo sample of it: f
Pj =
Nj
N , (2.35)
with Nj the number of particles that were absorbed and N the total number of particles in my simulation run.8 Thereafter, I can do statistical analysis on these data points.
8I usee here to denote the Monte Carlo sample, instead of the usual overbar for Monte Carlo
samples, because I am going to take the mean of my Monte Carlo samples later on, which is also denoted by an overbar.
2.3. Statistical methods 19 ABSORBED ESCAPED OUTSIDE CELL INSIDE SURFACE START d∗ RS LAYER
Figure 2.5.: A flow chart of the full algorithm with next to it all the parts of the simulation described in figures2.1to2.4.
2.3 Statistical methods
First, I want to determine the statistical error of my individual samples. That is, if I have some sample for the probabilistic influx fPj i, I would like to have an idea for how much the difference is with the actual probabilistic influx:
ei= Pj− fPj i. (2.36)
To make an estimate of this error,ebi, I need to determine the standard deviation eσ of the distribution of my samples fPj [21, pp. 1–2]. I note that my Monte Carlo
sim-ulation scheme is basically a Bernoulli process with the probability of the particle being absorbed Pj and the particle escaping the simulation 1− Pj. Therefore, the
number of particles that are absorbed, Nj, are binomially distributed. This is a
bi-nomial distribution with a meanhNji = NPj and a variance Var(Nj) = N Pj(1−Pj).
Assuming that Pj is not close to zero or one, and the number of particles in my
simulation is large enough; I can use the de Moivre–Laplace theorem [22], and proximate this binomial distribution with a Gaussian distribution. Hence, I can ap-proximate the distribution of my samples too with a Gaussian distribution, centred around Pj: P (fPj)' 1 √ 2πeσ2e −( fPj −Pj )2eσ2 2, (2.37) where eσ =pVar(Nj)/N =
p
Pj(1− Pj)/N is the estimate of the error I are looking
for.
I see that my estimate scales inversely with the square root of the number of particles in my simulation: eib = eσ ∝ 1/√N . I use a million particles in my simulations, and since the maximum value ofpPj(1− Pj) is1⁄2, I have thatebi≤1⁄2000. This means a
estimated by: σP j ≈ dσPj = s √ n. (2.39)
The standard error of the mean is estimated, because I do not know the true stand-ard deviation σ of the distribution I took samples from. Hence, I have to use the (corrected) sample standard deviation as an estimate for it:
s = v u u t 1 n− 1 n X i=1 f Pj i− Pj 2 , (2.40)
with 1/(n− 1) the Bessel’s correction. This correction decreases the bias in the estimate for σ, though it is still biased.
Using these statistical methods, I can now look at the data generated by my simula-tions.
3. Results
Before I could run simulations on different type of cell shapes with different distri-butions of transfer sites in their membrane, I had to check if my simulation method worked properly.
3.1 Validation algorithm
In order to validate my algorithm, I studied a case for which I do have an analytical solution. I took my absorbing body to be a sphere itself, with radius RC, and a homogeneous distribution of transfer sites, so that the absorption rate k(r) is some constant value k. Now equation (2.24) becomes:
Φj,cell = Φj,S2(0,R C)= Pabs S2(0, RC)|S2(0, RS) S2(0,R S)· Φjabs,S 2(0,R S), (3.1)
and I can use equation (F.13): Φj,S2(0,R C) = 4πRCc∞D kRC D + kRC , (3.2) and equation (F.14): Φjabs,S2(0,R S)= 4πRSc∞D, (3.3) so that: Pabs S2(0, RC)|S2(0, RS) S2(0,R S)→ Φj,S2(0,R C) Φj,S2(0,R S) =4πRC 2Dkc ∞ D + kRC 1 4πDRSc∞ =RC 2 RS k D + kRC , (3.4)
where the→ sign is used in the first expression to denote that the average of the probability that I will find has to converge to the exact number on the right-hand side, given that enough simulations are run.
I simulated ten sequences of a million particles each. Each sequence provided a sample (equation (2.35)), and I calculated the mean of those ten samples (equa-tion (2.38)). In figure3.1 on page22, I see that those means follow the analytical solution perfectly. The differences between the means and the analytical solution are shown explicitly in figures 3.2and 3.3on pages23 and 24, as the true error and the percentage relative true error, respectively. I see that for higher values of k there is a small but notable discrepancy. I have not been able to find the cause for this.
These disparities cannot be explained by statistical fluctuations. The estimates for the standard error of the mean (equation (2.39)), denoted by the error bars
0 2 4 6 8 10 Robin proportion k 0.0 0.1 0.2 Probabilistic influx Data simulation Exact solution (3.4)
Figure 3.1.: Comparison between the mean of the data generated by the algorithm and the exact solution for the probabilistic influx (3.4), as a function of the Robin proportion k, for a sphere with a radius RC = 2 and a uniform distribution of the Robin proportion k over its surface. The error bars for the data denote the estimate for the standard error of the mean, given by equation (2.39).
in figure 3.1on page 22, are minuscule. I present them explicitly in figure 3.4on page25. Here, I see that the error bars are in the order of 2·10−4. These diminutive error bars tell us that even though I do not have many observations — the mere ten sequences — my sampled means of the observations are still very close to their popu-lation mean. This closeness is caused by the high number of particles that are used per observation: one would expect each observation to converge to the same result due to the low random error. As the error bars are smaller than the true errors of my results (figure3.2 on page 23), the true errors are caused by the systematic error of the methodology and simulating more sequences would not improve them.
Regardless of the errors, figure 3.1on page 22 demonstrates that my algorithm is capable of capturing the correct behaviour of the system I am investigating. As such, the design is functional, and I trust that it will be so for other systems as well.
3.1. Validation algorithm 23 0 2 4 6 8 10 Robin proportion k 0.0000 0.0005 0.0010 0.0015 0.0020 T rue error figure 3.1
Figure 3.2.: The true error of figure3.1: the difference between the mean of the data generated by the algorithm and the exact solution for the probabilistic influx (3.4), as a function of the Robin proportion k, for a sphere with a radius RC = 2 and a uniform distribution of the Robin proportion k over its surface.
Positive numbers signify an overshoot of the data, negative numbers an undershoot.
0 2 4 6 8 10 Robin proportion k −0.1 0.0 0.1 0.2 0.3 0.4 0.5 P ercen tage relativ e true error figure 3.1 (%)
Figure 3.3.: The percentage relative true error of figure 3.1: the same difference (figure 3.2) between the mean of the data generated by the algorithm
and the exact solution for the probabilistic influx (3.4), now as a percentage of the exact solution. As a function of the Robin proportion k, for a sphere with a radius RC = 2 and a uniform distribution of the Robin proportion k over its surface.
Positive numbers signify an overshoot of the data, negative numbers an undershoot.
3.1. Validation algorithm 25 0 2 4 6 8 10 Robin proportion k 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 Error bar values figure 3.1
Figure 3.4.: Showing explicitly the order of magnitude of the values of the error bars in figure3.1. The values give the estimate for the standard error of the mean, given by equation (2.39), of the data generated by the algorithm for the probabilistic influx, as a function of the Robin proportion k; for a sphere with a radius RC = 2, and a uniform distribution of the Robin proportion k over its surface.
0 1 2 3 4 5 6 Thickness diffusion layer d∗(σ)
0.0 0.2 0.4 0.6 Probabilistic influx k = 1 k = 0.03
Figure 3.5.: The effect that the thickness of the diffusion layer d∗ has on the probabilistic influx of a sphere with radius RC = 1, for different values of k. The expected result from the exact solution (3.4) is shown in colour. Note that most data points shown here are the result of single simulations of a million particles, rather than the mean of ten of such simulations. Therefore, there are no error bars in this figure.
3.1.1 Thickness Gaussian layer
In order to find an optimal value for the thickness d∗ of the Gaussian layer, I con-ducted a parameter sweep in d∗ for different values of k, with the other paramet-ers fixed (I set ∆t to 0.0001 for precision in the Gaussian layer; the resolution is discussed in the end of section 2.2.2 A lattice random walk or a Gaussian random walk, on page 13). The results of this sweep are shown in figures 3.5and 3.6on pages 26 and27. Note that here the radii of the absorbing sphere, RC, and of the surrounding sphere, RS, differ from those of the previous figures.
I see that for very low values of d∗ in terms of σ, not enough particles are absorbed, as expected. The explication would be that it is difficult for particles to get out of the walk-on-spheres layer into the Gaussian layer. An unexpected finding is, that when d∗ increases, I start seeing overshoots. However, as described before, the al-gorithm does seem to grasp the correct behaviour of the system nonetheless. For my other simulations, I chose to set the thickness of the Gaussian layer to five times the standard deviation of a diffusive moving particle: d∗= 5σ.
3.1. Validation algorithm 27
0 1 2 3 4 5 6
Thickness diffusion layer d∗(σ) 0.80 0.85 0.90 0.95 1.00 Probabilistic influx normalized to exact solution k = 10 k = 1 k = 0.03
Figure 3.6.: The effect that the thickness of the diffusion layer d∗ has on the probabilistic influx of a sphere with radius RC = 1, now normalized to the exact solution (3.4) for different values of k. Note that most data points shown here are the result of single simulations of a million particles, rather than the mean of ten of such simulations. Therefore, there are no error bars in this figure.
ϵ 1− Pk
with D the diffusion coefficient, ϵ “is equal to the step size”, and Pk “represents the probability that a walker contacting the trap surface is absorbed”.
This absorption probability Pk is what I am interested in, so I can rewrite:
Pk= κϵ
D + κϵ
= κ
κ + (D/ϵ).
(3.6)
Attempting to interpret equation (3.6) for my situation, I set their κ to my k and put the step size ϵ to the mean-square-displacement of: √6D∆t. Therefore, I use:
Pk=
k
k + (D/√6D∆t). (3.7)
Vaughn proposed the following expression for the absorption probability [24]2: Pa= (1 + δ− x0/a)p
1 + δ(1 + p) x0
a, (3.8)
where a is the location of the absorbing surface (in one dimension, they consider a one-dimensional flux when the particle is about to be absorbed), x0 is the distance of the particle from a, δ is related to their diffusion layer, and p is the surface per-meability: p = κa/D, with κ the “rate constant of particle absorption” and D the diffusion coefficient.
When I try to translate equation (3.8) to my simulation, I again set their κ to my k; moreover, they state that aδ is “the thickness of the trapping zone”, hence for us: δ = d∗/RC. Doing this, I get:
Pa= (1 + δ− krik/RC)p 1 + δ(1 + p) krik RC , (3.9) with now p = kRC/D.
Comparing the absorption probability equation (2.27) with equations (3.7) and (3.9), resulted in figures3.7to3.9 on pages29 to31.
I see that Pabs(k) (2.27) yields better results for most values of k; it gets outper-formed by Pa (3.9) for k & 7, under these conditions — the radius of the absorbing sphere: RC = 1, of the surrounding sphere: RS = 1.1, the diffusion coefficient: D = 1, and the time increment of the Gaussian jumps ∆t = 0.0001. It is notable that Pa (3.9) generally performs quite well too. Regardless, these results strengthened
1Equation (7) there. 2Equation (11) there.
3.1. Validation algorithm 29 0 2 4 6 8 10 Robin proportion k 0.0 0.2 0.4 0.6 0.8 Probabilistic influx Exact solution (3.4) Data simulation Pabs(2.27)
Data simulation Pk(3.7)
Data simulation Pa(3.9)
Figure 3.7.: Comparison between the mean of the data generated with different expressions for the absorption probability: Pabs, Pk, and Pa
(equations (2.27), (3.7), and (3.9)), as a function of the Robin proportion k; for a sphere with a radius RC = 1 and a uniform distribution of the Robin proportion k over its surface.
The error bars for the data denote the estimate for the standard error of the mean, given by equation (2.39).
The exact solution for the probabilistic influx is shown as reference (3.4).
my confidence to use Pabs(k) (2.27) for my simulations of different absorbing bodies than the sphere.
0 2 4 6 8 10 Robin proportion k 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 T rue error figure 3.7
Data simulation Pabs(2.27)
Data simulation Pk(3.7)
Data simulation Pa(3.9)
Figure 3.8.: The true error of figure 3.7: the difference between the exact solution for the probabilistic influx (3.4) and the mean of the data generated with different expressions for the absorption probability: Pabs, Pk, and Pa (equations (2.27), (3.7), and (3.9)), as a function of the Robin proportion k; for a sphere with a radius RC = 1 and a uniform distribution of the Robin proportion k over its surface.
Positive numbers signify an overshoot of the data, negative numbers an undershoot.
3.1. Validation algorithm 31 0 2 4 6 8 10 Robin proportion k 0 5 10 15 20 25 30 35 P ercen tage relativ e true error figure 3.7 (%)
Data simulation Pabs(2.27)
Data simulation Pk(3.7)
Data simulation Pa(3.9)
Figure 3.9.: The percentage relative true error of figure3.7: the same difference (figure 3.8) between the exact solution for the probabilistic influx (3.4)
and the mean of the data generated with different expressions for the absorption probability (Pabs, Pk, and Pa; equations (2.27), (3.7), and (3.9)), now as a percentage of the exact solution. This difference as a function of the Robin proportion k; for a sphere with a radius RC = 1 and a uniform distribution of the Robin proportion k over its surface. Positive numbers signify an overshoot of the data, negative numbers an undershoot.
