Simulating Stochastic Reaction-Diffusion Systems on and within Moving Boundaries

Chemical reactions inside cells are generally considered to happen within fixed-size compartments. However, cells and their compartments are highly dynamic. Thus, such stringent geometrical assumptions may not reflect biophysical reality, and can highly bias conclusions from simulation studies. In this work, we present an intuitive algorithm for particle-based diffusion in and on moving boundaries, for both point particles and spherical particles. We first benchmark our proposed stochastic method against solutions of partial differential equations in appropriate scenarios, and further demonstrate that moving boundaries can give rise to super-diffusive motion as well as time-inhomogeneous reaction rates. Finally, we conduct a numerical experiment representing photobleaching of diffusing fluorescent proteins in dividing Saccharomyces cerevisiae cells to demonstrate that moving boundaries might cause important effects neglected in previously published studies of cell compartmentalization.


Introduction
Moving boundaries can have significant effects in the simulation of biological systems. For instance, it has been suggested that rapidly changing geometries in dividing cells might be important factors driving asymmetric cell division [1], while growing domains have already been identified as important elements in accurate models of cell migration [2] and tumour growth [3], among many examples.
Previous work in the area has mainly focussed on deterministic approaches, which might prove unsuitable for many biological applications. This is especially the case when the population of reactants can be on the order of only a few individuals, or when noise may propagate within the system causing non-classic effects [4]. In a stochastic setting, there are some examples of growing domains in the context of the reaction-diffusion master (Eq 2), where particles diffuse on a discrete lattice, and each particle is assumed to be evenly-distributed within each voxel on the lattice. However, there are significant drawbacks to the master equation approach: in the limit of a very fine lattice, master equations cease to account for bimolecular reactions [5]; and the discrete nature of the lattice fails to account for particle size, making it problematic to account for volume exclusion effects. Perhaps even more importantly for our purposes is that master equation approaches can introduce inaccuracies at boundaries [6] for several types of boundary conditions. Thus, there is an outstanding need for a simulation methodology incorporating both moving boundaries and stochastic effects, without the introduction of the artefacts that the reaction-diffusion master equation suffers from.
Alternative approaches for simulating stochastic reaction-diffusion systems with moving boundaries might draw from existing methodologies available for static domains. Potentially promising approaches which do not suffer from the aforementioned drawbacks treat space continuously, and can broadly be placed into two categories: simulators with fixed time steps, and simulators with adaptive time steps. For simulation techniques with adaptive time steps [7,8], the general approach is to consider imaginary domains around each diffusing particle such that the imaginary domain of one particle does not overlap with that of another. Analytical forms for the Green's functions for single particles can then be used to calculate the first hitting times to each imaginary domain boundary, and to propagate the particles at each time step chosen according to the smallest first hitting time. This allows the simulator to consider particles independently of each other at different time steps, in which case analytical solutions from single particle systems can be used. For interactions with boundaries, analytical forms for the Green's function with the appropriate boundary are needed. Extending such approaches to moving boundaries can become problematic, since Green's functions must be calculated for problems with moving boundaries, every time step. In general, these are not equations that lend themselves to analytical solutions, so numerical approaches become necessary for most simulation domains. Schemes for handling these problems do exist, such as moving finite elements [9], or Level-Set methods [10]. These approaches are commonly used to investigate problems involving solving partial differential equations (PDEs) across phase-changes, and they can be generalized to time-varying domains. However, these numerical approaches can entail significant computational costs, which makes using them unfeasible, especially when simulating systems with many particles. It is possible to choose an adaptive time step such that at most one particle interacts with a boundary per adaptive time step, thus limiting the required number of PDEs with moving boundaries to be calculated to a single instance. However, this would limit the adaptive time step to take small increments, and by consequence greatly hinder the numerical efficiency of the algorithm.
In this paper, we propose that the most promising approach for simulating reaction-diffusion in and on moving boundaries in feasible computational timeframes is to build on the basis of particle-tracking methods [11]. The crux of the algorithm is to choose a fixed time step, and propagate particles according to distributions consistent with their Green's functions in freespace. If a particle should fall outside the domain after a time step, it can be mapped back into the domain according to some rules. While such approaches have been used in static domains [11], we address here their applicability to dynamically changing domains.
We demonstrate the accuracy of our approach by numerically showing that particle-trackers can give sample paths consistent with corresponding PDEs in scenarios with moving boundaries. This investigation is conducted for point particles, as well as spherical particles with finite radii. Additionally, we investigate the effect of moving boundaries on mean-squared displacements, and how they can give rise to unexpected, non-standard reaction rate profiles. As a final application, we numerically investigate photobleaching in a dividing S. cerevisiae cell. This popular experiment has so far only been investigated in static domains, and we demonstrate here that moving boundaries can induce important effects that had not been accounted for so far.

Formulation
We assume the probability of a particle to simultaneously react and collide with a moving boundary to be negligibly small, thus allowing us to consider particle-boundary and particle-particle events separately. With this in mind, our aim is to construct a simulator that generates sample paths consistent with a suitable partial differential equation (PDE), i.e. the diffusion equation, for the probability density of a single particle with moving boundaries. After constructing this diffusive motion, interactions between particles may be added by established techniques.
For the following sections, we consider a closed time-dependent domain O t 2 R d (d 2 N) with a smooth boundary Γ t . We denote the region outside of the boundary by O c t . For every x 2 Γ t , we denote the velocity and the outward normal vector to the boundary by v t ðxÞ ¼ ðv ðtÞ 1 ; v ðtÞ 2 ; . . .; v ðtÞ d Þ and n t ðxÞ ¼ ðn ðtÞ 1 ; n ðtÞ 2 ; . . .; n ðtÞ d Þ, respectively. Note that we consider a ddimensional Euclidean space so as to provide a general setting for our algorithm.
Diffusion and reaction events within moving boundaries. We assume that the probability density function (PDF), P, of a single particle in O t undergoing Brownian motion obeys the following relation: where the boundary condition has been derived by imposing the conservation of probability. Namely, assuming that the boundary condition takes the form rP = g(P, v t , n t ) for all x 2 Γ t , the following must hold: where the first equality follows by applying Leibniz's integral rule. The second equality follows directly from the Divergence Theorem, namely by equating the rate of change of probability in the domain with the flux of probability through the domain boundary. Comparing terms yields the boundary condition in Eq 1. It is worth noting that for a static domain, v t (x) = 0 for all x 2 Γ t holds, and the boundary condition takes the form of a zero-Neumann type. In what follows, we will concentrate on Neumann boundary conditions. However, generalisations of the particle tracker to incorporate effects of Robin boundary-conditions directly follow from existing methods of semi-permeable membranes in the case of static boundaries [11]. Since Eq 1 represents a diffusion equation with particle conservation within the domain, a natural simulation procedure would be to propagate each particle according to a Gaussian distribution, and map any escaped particles back into the domain. Thus we constructed an algorithm represented by the pseudo-code in Fig 1 and graphically in Fig 2a. In summary, we introduce n 2 N particles, which we enumerate with some index i 2 {0, 1,. . ., n}, each with radius ρ i and associated diffusion coefficient D i . Note that the explicit incorporation of particle size allows this method to circumvent problems associated with master equation (please refer to the Introduction for further details).
Each particle is propagated in a small time step, dt, along each Cartesian coordinate according to a Normal distribution with mean zero and variance 2D i dt. At each time step t, we check for any particles lying outside O t or within the particle radius away from Γ t . For each of these particles, we calculate the point on Γ t which is closest to the escaped particle, and construct a vector, r, from the particles position, x i , to this point. For point particles (i.e. ρ i = 0), we set x ðnewÞ i ¼ x i À kr i for some k > 1. For particles with positive radius, we place the particles at the closest point on the domain, and then move them in the normal direction by ρ i . This ensures that the particle is contained within the domain.
Regarding chemical reactions, we execute a loop over all pairs of particles to calculate the distance between each possible pair of particles. If reactive particles i and j are separated by a distance less than or equal to l ðrÞ ij , we execute reaction r accordingly. In case there is a single reaction product, we choose the location of the product to be at the same position (unary reactions) or in the middle position between the two reactants (binary reactions). In case there are two products, we place the particles in the same locations as the original reactants for the case of binary reactions, or at randomly opposite poles of a small sphere centred on the reactant in the case of unary reactions. The size of the sphere is chosen to be of the order of the largest reaction radius,l ðrÞ ij . In some scenarios, this might cause artefacts (for instance, in case of 'germinate' reactions, where reaction products react with each other immediately). Further artefacts might occur since we do not consider particles reacting between time steps. In these cases, and in cases for a higher number of reaction products, alternative procedures from static boundary scenarios can be used [11]. For the case that two particles' radii overlap with each other after a time step, the particles are then moved an equal amount away from each other in the direction defined by the vector connecting the particles' centres. Note that the algorithm represented above allows for a variety of effects at the boundary. For example, setting k = 1 amounts each escaped particle to be placed at the closest point on the domain boundary, and using k = 2 amounts to reflecting a point particle through the closest point on the domain boundary to a point within the domain boundary. We seek to demonstrate first that, for point particles without reactions, the sample paths generated from our algorithm are consistent with Eq 1 for a variety of different boundary effects. Please see subsequent sections for benchmarking of the algorithm.
Diffusion and reactions on a moving boundary. We assume that the PDF, P, of the position for a diffusing species on Γ t obeys the following equation: where r 2 G t ¼ ðr À n t ðn t :rÞÞ:ðr À n t ðn t :rÞÞ refers to the surface Laplacian. Note that if P is initially defined on Γ 0 , then the advection term in Eq 3 will act in a way such that P remains on Γ t for all t. Also note that Eq 3 represents diffusion on a surface with an advective term in the direction of the surface movement. Thus, a logical approach to generate sample paths consistent with Eq 3 is to propagate particles according to a Gaussian distribution at every time step to account for motion tangential to the surface, and further add in an advective term given by the surface motion. Accordingly, we construct the algorithm represented by the pseudo-code in Fig  3, and graphically in Fig 2b. As before, we have n 2 N particles with the same index set i 2 {0,1,. . .,n}, with particle radii ρ i and diffusion coefficients D i . Particles i and j react in a similar way as the previous section, i.e. when their centres are separated by a distance of less than l ðrÞ ij . At every time step, each particle is propagated in a small time step, dt, along each Cartesian coordinate according to a Normal distribution with mean zero and variance 2D i dt. For every time step from t to t+dt, the particles are further propagated by the amount v t, dt. For each particle, we then calculate the point on Γ t which is closest to the escaped particle, and map each particle onto this point. Then, we execute reactions and collisions exactly as per the case for diffusion within the volume. Note that this procedure is independent of the curvature of the surface of the boundary. A justification for this follows from considering the limit of small time steps, such that the diffusive motion is small, and the surface is locally flat. Since the projection of a Brownian motion itself is Brownian, this procedure should give consistent sample paths to Eq 3.

Algorithm Benchmarking and Sample Theoretical Applications
Diffusion within and on a shrinking circle. We first consider diffusion within and on a shrinking circle as a sample application. We first demonstrate that the sampled particle tracks generated with our algorithm are consistent with the solutions of the PDEs in Eqs Eq 1 and Eq 3, then proceed to investigate the effects of moving boundaries on the mean squared displacement (MSD) of diffusing particles. Note that we treat PDEs as descriptions of the PDFs of individual particles, as opposed to the behaviour of diffusion of many particles.
Consider O t to be the area bound by a circle in R 2 , centred at the origin, with Γ t denoting its boundary. At time t, we denote the radius of such a circle by r t . The starting radius is written r s , and the finishing radius as r e . Simulations occur up to an end time of t e , and the change of r t from the starting radius to the end is taken to be linear in time. We denote the Euclidean distance of a particle at time t from its initial position by x t . For hx 2 t i / t a , we denote the process as anomalous diffusion: of the sub-diffusive type for α < 1 and super-diffusive for α > 1. Diffusion is denoted as Brownian, or non-anomalous, if α = 1. In all subsequent analyses, to estimate the error in the empirical CDFs from simulated sample paths, we make use of a normal approximation to the binomial distribution, and employ a confidence interval of 99%. That is, for some empirically estimated, binomial probabilityp, we estimate the error as AE2:56 , where n is the number of simulated trials. The factor of 2.56 is established by considering the 99.5 percentile of a standard normal distribution.
To estimate the error in MSDs, we make use of the standard error of the mean. That is, for an empirically determined standard deviation, s, of a population of n sampled squared distances, we estimate the error in the mean as s ffiffi n p . Diffusion inside the volume of a shrinking circle. We numerically demonstrate convergence by considering the cumulative density functions (CDFs), P(R < r t ) of the radial component (R) of diffusing particle. We consider shrinking circle domains with uniform initial conditions, incorporating diffusion constants in widely different scales. Fig 4 shows results considering several scenarios: 1) the diffusive motion of each particle is much faster than that of the boundary (Fig 4a and 4d), 2) particles and the boundary moving at a similar rates (Fig 4b  and 4e) and, 3) particles diffusing much slower than the boundary (Fig 4c and 4f). We investigate these three scenarios for both point particles (Fig 4a-4c) and for particles with finite size (Fig 4d-4f). For point particles, we consider two methods by which to map escaped particles back into O t : we place escaped particles at the closest point on the boundary (i.e. k = 1 in the pseudo-code in Fig 1), or we reflect the particles through the closest point in the boundary (i.e. k = 2). The error in CDFs between the two methods were found to decrease as simulation time steps were reduced, as demonstrated in Fig 4g. Solutions to Eq 1 were found adopting a deforming mesh from the commercial software package Comsol Multiphysics version 4.3. To allow for particle size to be included into the PDE solutions, the value r t was reduced by an amount ρ for all t. Stochastic simulations were run with single particles with initial uniformly distributed positions, and were repeated 10000 times. Initially, simulations were repeated 2000 times, and indistinguishable results to those with 10000 repetitions were obtained. Thus, the number of repetitions was deemed sufficiently large to derive conclusions for this part of the study. Please refer to the legend of Fig 4 for the time step of stochastic simulations. In these scenarios, we find that the CDFs generated by the particle tracker and PDEs are in good agreement, and this is independent of the correction method used to map escaped particles.
In cases where the diffusion is much faster than the boundary movement, we observe a uniform distribution in the circle at all times (i.e. P(R < r t ) = R 2 / (r t −ρ) 2 ). In situations where diffusive motion is significantly slower than the boundary movement we observe that the CDFs of particles remain consistent with the uniform distribution over the initial circle with radius r s , rising quickly to unity close to the boundary of the sphere i.e.: ( These analytical results are in good agreement with numerical experiments, as shown in Fig  4c and 4f. This model system also allows us to consider the effects of moving boundaries on the mean squared displacement of diffusing particles. For doing this, we denote the Euclidean distance of a particle at time t from its initial position by x t . For large diffusion coefficients, we observe that the mean-squared displacement is consistent with freely diffusing particles (hx 2 t i / t) so long as we consider short time scales to reduce the effects of diffusion within confined volumes  (note the coefficient close to unity of the linear regression represented by red dots in Fig 5c). At larger times, the beginning of confinement effects give rise to sub-diffusive behaviour [12] (note the black dots in Fig 5c). Assuming the particles are diffusing very slowly relative to the boundary, we conjecture that the MSD can be decomposed into two components: the fraction of freely diffusing particles, and those being pushed by the boundary. Enumerating these terms results in the following equation, where the first term represents freely diffusing particles, and the second term represents those being swept from their initial positions by the boundary: In more detail: The first term represents the product of the proportion of freely diffusing particles ( Note that such arguments so far would neglect any effect of growing domains, i.e. where the boundary would be moving away from particles. In such scenarios, we would expect our current methodology to give non-anomalous MSD characteristics. However, it may well be that such boundaries give rise to an advective velocity on moving particles and might still provide interesting effects. In such scenarios, advective forces can be introduced into particle tracking methodologies, but we consider this is outside the scope of the current work. Diffusion on the surface of a shrinking circle. We now consider the case of diffusion on the surface of a shrinking circle, with initial conditions as a point mass at the coordinate with the largest y-value in an x-y Cartesian plane. In this scenario, numerical solutions of Eq 3 can be found by considering the problem in polar coordinates. We perform the transformation in two steps, first by considering a transformation to arc-lengths along the circle (which removes advective terms due to the moving boundary), then to arc-angles (which transforms the problem from a moving boundary to a static boundary problem). We first define the position of a particle in terms ofr t 2 ðÀpr t ; pr t , which is the arc-distance along the shrinking circle of the particle from the point where the circle intersects the positive y-axis. We taker t to be positive (negative) for coordinates with a positive (negative) x coordinate, to denote the specific displacement. No advective term results in the transformed diffusion equation, since the movement of the surface is perpendicular to the diffusive motion along surface. Now we can reformulate the problem in terms of a diffusion equation for the PDF Pðr t ; tÞ of the particle being atr t at time t: The above equation is a moving boundary problem, but can be reformulated into a static domain by a change of variables to the angle from the positive y axis, here denoted by θ 2 (−π, π], and we write the corresponding PDF asPðy; tÞ. Using the relationshipr t ¼ r t y, we can transform Eq 6 into the following equation. which is a diffusion equation with a time-varying diffusion coefficient, but with periodic nonmoving boundaries. This formulation lends itself to a simple numerical solution via an explicit finite different method, and allows us to benchmark the proposed algorithm. By considering N equally spaced points θ 1 ,θ 2 ,. . ., θ N spanning [−π,π), with grid spacing of h ¼ 2p NÀ1 and time steps of δt, we construct the PDE solution in terms of the probability density at point θ i , i 2 {1,. . ., N} at time t, which we write as P t y i . For all i 2 {1,. . ., N}, we use the following scheme to generate PDE solutions: where, if i = 1 we define i-1 N and if i = N, we define i + 1 1to signify that the domain is circular.
Numerical investigations over a range of diffusion coefficients show that our algorithm generates consistent sample paths with Eq 7, as shown in Fig 6a-6c. We further investigated the effect of the moving boundary on the MSD. At large diffusion coefficients, the simulated behaviour is consistent with free diffusion at the start, before becoming sub-diffusive as the MSD approaches an upper bound (due to the finite circumference of the circle at the end of the process). To derive an expression for the upper bound of the MSD, assume that the distribution of particles is uniform over the circle at the end time. Considering the initial particle position of (0, r s ), and a final position (r e cosθ, r e sinθ) that is uniformly distributed over θ 2 [0,2π), this gives the following expression: ð 2p 0 ððr s À r e cosyÞ 2 þ r 2 e sin 2 yÞdy ¼ r 2 s þ r 2 e : ð9Þ In the limit of zero diffusion rates, we can derive the limit behaviour by assuming the diffusing particle is static on the moving domain at the point (0,r t ) for all t. This trivially gives the expression: Recalling that r t changes linearly in time, this implies that the motion is super-diffusive if the displacement due to the moving boundary is larger than that related to particle diffusion. As diffusion rates increase such that the motion due to diffusion is of comparable to that of boundary movement, we observe that super-diffusive effects subside, as shown in Fig 6d. Now that we have presented our algorithm in full and benchmarked it against known solutions, we will focus on sample applications to illustrate its ease of use and applicability.

Sample Theoretical Applications
Bimolecular decay inside an elongating dumbbell. We now consider some non-standard effects that moving boundaries might have on reaction kinetics. The domain in consideration is two spheres joined together by a bridge, representing an elongating dumbbell. This shape resembles dividing nuclei in closed mitosis [13], among other biophysical processes and structures. Initially, the bridge is of length zero. After some time, the length of the bridge increases steadily toward a final value, upon which it stops growing. A schematic of the final shape of the sphere projected in 2D is found in Fig 7. Thus, we have a domain in R 3 with rotational symmetry around the x-axis. We consider a particle to be within the domain so long as the position of each particle (x(t), y(t), z(t)) satisfies: . This system represents two spheres of radius R, joined together by a cylinder of diameter d with length l(t), a geometry common in dividing cells/ nuclei. The centre of the cylinder is at x = 0.
We then consider bimolecular decay A + B ! ;, where particles of species A and B react when they are within a distance λ at the end of a time step. In this example, the probability of reacting is always 1 whenever particles are within that predefined distance. However, the algorithm can be generally set up in such a way that the reaction probabilities take values in the interval [0,1] once a distance threshold is reached. As an initial condition, we place a single A molecule at the centre of the left sphere, and a single B molecule at the centre of the right sphere. We then define the length function l(t) as follows: i.e. the length is monotonically non-decreasing. In a fixed volume, one expects the reaction times to follow an exponential distribution. However, in our moving domain scenario, even though the total volume of the system is increasing, we note that the time-distribution of the decay event can be bimodal, as shown in Fig 7. This interestingly implies the effective reaction rate can increase with increasing volumes. In Fig 7, we allow the bridge to remain at length zero for 200,000 seconds, before allowing the bridge to grow uniformly to reach length 10 units at 201,000 seconds. The spheres have a radius of 3 units at all times. Each particle diffuses slowly, with a rate of 0.0025 units squared per second. Particles react when they reach a distance of less than 0.4 units of each other. The resulting histogram was constructed from 50,000 events, with the bin-size determined by the Freedman-Diaconis rule [14]. As the bridge elongates, we observe an increased frequency of reaction events inside this region. Simulations with increasingly small time steps ensure that the bimodality in the histogram is not due to numerical error. For illustration purposes, Fig 7 shows a superposition of the geometry (projected onto the x-y plane) at 201,000 seconds on all reaction events occurring between 200,000 and 201,000 seconds. Note that only the final geometry is depicted in Fig 7. By consequence, some reaction events occurring before the geometry reaches its final conformation appear outside the domain.
Photobleaching in and on a dividing budding yeast nucleus. As a final application, we consider photobleaching experiments done in a dividing S. cerevisiae nucleus during anaphase. Note that S. cerevisiae undergoes closed mitosis, whereby the nucleus remains intact throughout anaphase, only to break down moments before cytokinesis. Our numerical simulations aim to reproduce a more physically realistic scenario of S. cerevisiae nuclear division than those considered in previously studies [13,15].
We represent mother and daughter nuclear lobes by prolate ellipsoids connected by a cylindrical bridge, as depicted in Fig 8. The major axes of the mother and daughter, as well as the axis of the cylinder, are taken to be co-axial to the x-axis. Initially, we construct the dividing nucleus in Cartesian coordinates (x', y', z') 2 R 3 , with the extreme edge of the mother lobe coinciding with x' = 0, followed by the bridge, followed by the daughter lobe. Then we transform the coordinate system such that the centre of mass of the nucleus lies in the centre of a coordinate system (x, y, z) 2 R 3 . We denote the major and minor axes of the mother and daughter by 2r xm,t , 2r ym,t and 2r xd,t , 2r yd,t respectively. We take the ellipsoid to be rotationally symmetric around the x-axis. Furthermore, we define the length of the bridge to be l t and its width (i.e. radius) to be w t .
Separately, images of dividing nuclei were obtained by confocal microscopy on a Zeiss LSM-780 microscope, courtesy of Dr. Jai Denton. The strain visualized, GA-1320 containing the pAA4+lacO plasmid, was kindly provided by Professor Susan Gasser and has been previously described in [1]. For imaging, live S. cerevisiae cells were grown overnight in 5ml SC-hisura medium, diluted 400ul in 5ml of fresh SC-his-ura medium and grown for a further 4 hours before being harvested and immobilised between a 3% low melting temperature agarose gelpad and a coverslip for visualisation. The confocal images and the idealised geometry constructed from them can be found in Fig 9. Geometrical parameters were measured manually with Image J from these images, with time intervals of 50 seconds between successive frames. For intermediate times between frames, values were linearly interpolated to allow the shape to change smoothly. Each ellipsoid was truncated such that the aperture of the bridge and the ellipsoid coincide without any gaps. Thus, particles in the mother cell have an x-component that satisfies 0 x 0 < r xm;t þ r xm;t cos arcsin w t r ym;t , with the arcsin being taken at the smallest positive value. The set of all points in the mother lobe can be then described as: ðx 0 ; y 0 ; z 0 Þ ¼ ðr x;t cosy þ r xm;t ; r y;t sinysin; r y;t sinycosÞ; ð13Þ with r x,t 2 [0,r xm,t ], r y,t 2 [0,r ym,t ], ϕ 2 [0,2π) and y 2 arcsin w t r ym;t ; p h i . Note that the limits for x' in the mother lobe stop short of 2r ym,t , since a small amount of the ellipsoid is truncated so as to attach the bridge. The area in the bridge consists of those points for which the following relationship holds r xm;t þ r xm;t cos arcsin w t r ym;t x 0 < r xm;t þ r xm;t cos arcsin w t r ym;t þ l t : ðx 0 ; y 0 ; z 0 Þ ¼ ðx 0 ; w t sin; w t cosÞ ð 14Þ with ϕ 2 [0,2π). Finally, the points in the daughter lobe are those points with x' coordinate satisfying both x 0 ! r xm;t þ r xm;t cos arcsin w t r ym;t þ l t , and x 0 < r xm;t þ r xm;t cos arcsin w t r ym;t þ l t þ r xd;t þ r xd;t cos arcsin w t r yd;t where the largest value less than π is taken from the arcsin and: ðx 0 ; y 0 ; z 0 Þ ¼ r xm;t þ r xm;t cos arcsin w t r ym;t ! ! þ l t þ r xd;tþ r x;t cosy; r y;t sinysin; r y;t sinycos !
This geometry is then translated parallel to the x-axis such that the centre of mass is situated at the origin for all times. The geometry is updated at every time step of the simulation, ensuring a smoothly changing geometry and thus minimising numerical artefacts that might arise from more sudden changes of geometry. The velocity field of the moving boundary is calculated numerically, and is assumed to act in a way such that the velocity carries each point on a Fig 8. A schematic of the idealisation of a dividing yeast nucleus. Parameters r ym,t and r xm,t represent the maximal radius of the mother nuclear lobe in the direction of the y and x axis at time t, respectively. We define r yd,t and r xd,t analogously for the bud nuclear lobe. The length and width of the bridge is given by l t and w t , respectively, and the radius of the spherical photobleaching spot is given by r photo . Note that the geometry is in 3D, and only a 2D projection is depicted here.
doi:10.1371/journal.pone.0133401.g008 membrane to its nearest point at the subsequent time. This assumption might not necessarily hold true for dividing S. cerevisiae nuclei, but such data is difficult to measure empirically, and so here we must resort to some assumption of this nature.
Finally, we consider a 'bleaching' region in which a degradation reaction A! k deg ; occurs, where k deg gives the rate per unit time that a single protein A is photobleached. The region where photobleaching can occur is taken to be the volume of a sphere of radius r photo centred at (r xm,t , r ym , t , 0). That is, on the edge of the centre of the mother lobe, and such region will be denoted as the 'bleach spot' to conform with terminology used in [13,15]. At the end of each time step of length dt, any particles lying inside the photobleaching region can decay according to a Bernoulli trial with probability k deg dt.
To compare the effect of moving boundaries, we conducted two sets of simulations: one set with static boundaries, and another where boundaries were allowed to change with time. These sets of simulations were repeated for proteins diffusing both in the membrane and on the membrane. See S1 and S2 Videos for diffusion within and on moving boundaries, and S3 and S4 Videos for diffusion within and on static boundaries. We define the middle of the bridge to mark a dividing plane where proteins switch between mother and daughter lobes. Plots for the number of A proteins in the mother and daughter lobes are shown in Fig 9. These are split between diffusion on the nuclear membrane (Fig 9a) and diffusion inside the nucleus (Fig 9b). Three runs of simulations with 300 proteins per run were conducted, with proteins initially uniformly distributed through the whole nucleus. In the case of diffusion on the nuclear membrane, we find that the population of A proteins in the mother lobe remains similar between the static and growing cases. However, the growth of the domain allows a significant proportion of proteins to avoid photobleaching by diffusing/escaping to the daughter lobe. We observe that the average numbers of unphotobleached proteins on the entire nuclear membrane are 80 and 122 for the case of static and moving boundaries, respectively. Thus, the effective photobleaching rate for diffusing particles in this static scenario is higher than that of the moving boundaries scenario. Thus, neglecting this effect might lead to underestimates of bleaching rates in FLIP and FRAP studies. For photobleaching experiments performed for particles diffusing inside the nucleus, we observe a similar overall photobleaching rate between the static and moving boundaries cases. However, the distribution of proteins is markedly different: decay in the mother lobe is significantly faster and the growth of the domain causes a significant proportion of proteins to also move into the daughter lobe. Thus, we expect that moving boundary effects might have important influence over previously published estimates of compartmentalisation effects in living cells, which highly depend on the ratio of fluorescence decay rates between the mother and daughter lobes [13,15,16].

Discussion
We have numerically demonstrated that intuitive correction for moving boundaries in particle-tracking based simulators can give rise to sample paths consistent with PDEs representing the diffusion equation with moving boundaries, and that such effects can lead to non-standard effects such as super-diffusion and the creation of unusual reaction profiles. We considered both cases of particles diffusing within the domain boundary, and on the boundary. While there have been previous studies showing that master equation approaches can incorporate growing domains [2], many drawbacks of master equations, such as artefacts in bimolecular reactions and boundaries, have been well-documented [5,6]. These drawbacks motivated the approach we present here, based on particle tracking. Alternative approaches might look to extend the work of Green's Function Reaction Dynamics or First Passage Monte Carlo kinetics, but such techniques would involve numerically solving multiple moving boundary PDEs which, whenever feasible, would at least cause a significant computational overhead when compared to our approach presented here.
So far, we have chosen to disregard potential advective forces arising from boundary movement. Given the dynamic nature of biological systems, with active and passive transport across membranes, it is not possible to arrive at suitable analytical forms for advective forces arising from moving boundaries in a general setting. However, if such forces could be experimentally determined, then our presented approach could generalise to incorporate them.
The particle-tracking algorithms for moving boundaries were benchmarked against solutions of corresponding PDEs for the case of a shrinking circle, where we considered scenarios ranging from those where the diffusive motion was slow compared to boundary movement, to those where the diffusive motion is significantly faster in comparison. For diffusion within the volume, we find that the probability distribution of quickly diffusing particles can be well approximated by a uniform distribution at all times. The effect of the moving boundary is most noticeable at smaller diffusion coefficients, and thus we suggest that some of the most fruitful applications of moving boundary algorithms would be cases where particles diffuse slowly relative to moving boundaries. A numerical investigation into the MSD of moving particles in the static domain found that fast diffusing particles retained a profile resembling diffusion in a static domain. However, slower diffusing particles were pushed by the boundary, leading to scenarios with motion described by super-diffusion. Similar results stem from the case of diffusion on the moving boundary, where slow diffusing particles are found to show super-diffusive behaviour.
By considering the limiting case of very slow diffusion, we were able to derive analytical forms for the MSD curves corresponding to diffusion on and in the shrinking circle, and found these results to be consistent with our numerical results. We anticipate that such results should carry into higher dimensions as well, with some differences. As the integral in Eq 5 shows a dependence on dimensionality, we anticipate that the degree of anomalousness of diffusion should depend on dimension. These results might motivate alternative explanations for currently documented cases of superdiffusion [17], where our observations of superdiffusion on a moving membrane might seem pertinent.
The case of bimolecular decay in an elongating dumbbell was used to demonstrate some potentially unusual effects that moving boundaries might produce in non-linear systems. In particular, even though the overall volume of the dumbbell is increasing with time, we observe a bimodal distribution of reaction times. This implies that the reaction rate increases at some point in time. The reason for this is that, as the dumbbell elongates, the concentration of reactants within the bridge of the dumbbell increases and, consequently, an increase in reaction rate follows. Thus, there might be reason to expect novel effects of moving boundaries in other scenarios where non-linear effects manifest themselves. Indeed, some work has already been done on investigating Turing patterns on growing domains [18], and this remains a fruitful ground for further work. Furthermore, we expect that the movement of any boundary in the direction of reactants will cause an increase in the reactant concentrations near the boundary, thus leading to an increase in reaction rates.
Finally, we concluded our sample applications by investigating photobleaching experiments in dividing S. cerevisiae nuclei. Such experiments have been numerically conducted on the assumption of static domains, but we demonstrate that moving boundaries can have important effects which can give rise to noticeably different results. In particular, we observed that the effective photobleaching rate for particles diffusing on nuclear membranes is higher for static boundaries than for moving boundaries, implying that existing methodologies might be underestimating photobleaching rates. Of equal importance, in formulating the model of a dividing nucleus, we were forced into making assumptions on the nature of the velocity field induced by a moving nuclear membrane. In the absence of any empirical studies to guide us in the choice of such a velocity field, we constructed a velocity field on the basis of mapping points from the domain at one time to the closest points on the domain at another time. While mathematically convenient, it is not clear whether budding yeast nuclear membranes specifically move in such a way, and this motivates further experiments into the movement of nuclear membranes during anaphase. It would be interesting to conduct further empirical research such that any and all moving boundary effects could be fully incorporated into future models. However, the current study strongly suggests that changing geometries should be nevertheless included into photobleaching studies.
In conclusion, we find that intuitive based corrections to particle trackers give rise to sample paths consistent with the diffusion equation with moving boundaries. The effects of moving boundaries are diverse, ranging from causing super-diffusion in slow diffusing particles, to manipulation of local concentrations of reactants. Up to now, there has been relatively little work done in stochastic reaction-diffusion systems enclosed by moving boundaries, and it remains a fruitful ground for further research.