Multiscale modeling of diffusion in a crowded environment

Molecules undergoing diffusion and reactions inside living cells are often modeled by the reaction-diffusion equations. These are continuous, deterministic partial differential equations (PDEs) describing the time evolution of the concentrations of molecules. For a diffusing molecule the concentration $u(\mathbf{x},t)$ is described by the diffusion equation

$$u_{t}(\mathbf{x},t)=\gamma_{0}\Delta u(\mathbf{x},t),\quad\mathbf{x}\in\Omega$$

(2.1)

with diffusion coefficient $\gamma_{0}$ and suitable boundary conditions on $\partial\Omega$. To include reactions corresponding terms are added. This macroscopic description is accurate in the limit of large molecule numbers, when stochastic fluctuations are small and the mean value is the quantity of interest. Important molecules such as DNA or transcription factors are, however, only present at very low copy numbers inside living cells. It has been observed in experiments [18, 51, 56, 58, 64, 72] and shown theoretically [24, 52] that stochastic fluctuations play an important role and a discrete stochastic description is more accurate than the deterministic equations.

We distinguish two levels of accuracy of stochastic models. In the mesoscopic model the domain $\Omega$ is partitioned into $N$ non-overlapping voxels $\mathcal{V}_{i}$ with nodes $\mathbf{x}_{i}$ at the center. The state vector $\mathbf{y}(t)$ contains the number of molecules $y_{i}(t)$ in each voxel $\mathcal{V}_{i}$ at time $t$. The voxels are small enough that the molecules can be considered well mixed inside so that reactions can occur between molecules residing in the same voxel. An individual molecule can jump from a voxel $\mathcal{V}_{i}$ to a neighboring voxel $\mathcal{V}_{j}$ to model diffusion. The diffusion master equation (DME) describes the time evolution of the probability to be in state $\mathbf{y}$ of a system with only diffusion

$$\frac{\partial p(\mathbf{y},t)}{\partial t}=\sum_{i=1}^{N}\sum_{j=1}^{N}\lambda_{ij}(\mathbf{y}-\boldsymbol{\mu}_{ij})p(\mathbf{y}-\boldsymbol{\mu}_{ij},t)-\lambda_{ii}(\mathbf{y})p(\mathbf{y},t),$$

(2.2)

where $\lambda_{ij}$ is the jump propensity from $\mathcal{V}_{i}$ to $\mathcal{V}_{j}$ and $\lambda_{ii}=\sum_{j=1}^{N}\lambda_{ij}$ is the total propensity to leave voxel $\mathcal{V}_{i}$. The transition vector $\boldsymbol{\mu}_{ij}$ is zero except for $\boldsymbol{\mu}_{ij,i}=-1$ and $\boldsymbol{\mu}_{ij,j}=1$. Let $\theta_{ij}$ be the splitting probability, that a jump from $\mathcal{V}_{i}$ goes to $\mathcal{V}_{j}$, then

$$\lambda_{ij}=\theta_{ij}\lambda_{ii}.$$

(2.3)

By including reaction terms in a similar manner, the DME can be extended to the reaction-diffusion master equation (RDME). In the presence of bimolecular reactions there exists no analytical solution and a numerical solution is impossible due to the high dimension of $\mathbf{y}$. Instead, one samples trajectories of the system with the stochastic simulation algorithm (SSA), first presented by Gillespie [26] for only reactions and improved in [8, 25]. The algorithm was extended to space dependent problems with a Cartesian partioning of the domain in [15] implemented in [34]. The propensities $\lambda_{ii}$ are here used to generate random numbers for the time until the next jump. To represent the complicated geometries present in cells the algorithm was extended to curved boundaries in [40] and adapted for unstructured meshes in [19] with software in [14, 38].

In the more accurate microscopic model the molecules are tracked along their Brownian trajectories in a continuous space, continuous time Markov process. One approach is called Green’s function reaction dynamics (GFRD), methods and software for this approach are found in [1, 12, 43, 77] with a review in [69]. In [13, 61, 73] protective domains are constructed around individual molecules in which they cannot interact with other molecules. The exit times and exit positions from these domains are sampled to propagate the system until molecules are close enough to interact without sampling all the intermediate jumps.

The macroscopic, mesoscopic and microscopic models presented above are designed to simulate diffusion in a dilute medium. The microscopic model incorporates crowding effects automatically since the molecules are modeled as hard spheres with a given volume, but it becomes computationally very expensive in a densely packed space of inactive crowders, because the protective domains around molecules will be small and many short jumps will be simulated before meeting a potential reaction partner.

Cellular automata (CA) have been used in [5, 9, 68, 74] to simulate diffusion in a crowded environment. This is a lattice, or voxel, based microscopic approach, where each site can hold one molecule and crowders are represented as already occupied lattice points. The jump length is here the size of a molecule which also leads to expensive simulations with many short jumps and the shape of the molecules corresponds to the shape of the chosen lattice. The choice of the lattice, moreover, influences the excluded volume effect [30]. A Cartesian grid leads to a stronger crowding effect than a hexagonal mesh in 2D and both overestimate the effect of crowding on the reaction rates compared to a BD approach.

In [67] a mesoscopic approach is used where each voxel can hold more than one molecule. After distributing immobile crowders it is decided which voxels are accessible and which are full. The crowders can move in [21, 75] and the jump propensity to an adjacent node is scaled by the number of available spaces in the target voxel. Macroscopic nonlinear PDEs are then derived in [21] to model diffusion in the crowded cell envrionment and the results are validated by physical experiments in [20]. This approach is extended in [62] to derive nonlinear diffusion equations modeling more complicated interactions than steric repulsion between the molecules. Similarily, the averaged occupied volume in the whole domain is used in [45] to rescale the jump propensities. These approaches at most take the averaged occupied volume in the target voxel into account and neglect the microscopic positions, the shape of the molecules and the surrounding medium. Hence only an averaged behavior is observed and the MSD is linear like for normal diffusion but with a reduced diffusion constant:

$$\langle\mathbf{x}^{2}(t)\rangle=(1-\phi)\gamma_{0}t,$$

(2.4)

where $\phi$ is the occupied volume fraction.

In this paper we present a novel multiscale approach to simulate molecular crowding. We will use the microscopic information of the crowders’ positions to recalculate the jump propensities on an overlying mesoscopic mesh. Instead of simulating diffusion in the detailed environment, as in BD, we will have to solve many PDEs on small subdomains resolving the microscopic positions of the crowders. This is similar to homogenization techniques used e.g. for flow simulations in porous media, [7, 49]. The crowding molecules can have any shape and this approach is especially useful when the crowded environment is stationary or evolves on a much slower time scale than the diffusing molecule, so that the jump coefficients can be precomputed and used for a long simulation time. This is reasonable since the macromolecules responsible for the majority of occupied volume are ribosomes, microtubules and actin filaments [17], which are large and hence diffuse on a slower time scale than for example transcription factors. Moreover, it was shown in [11] that anomalous behavior is most likely to happen in a stationary environment, since otherwise averaging effects simply reduce the coefficient of normal diffusion. But, it is important to mention that by computing statistics our method can be made efficient also for moving crowders.

Let $c(\mathbf{x},t)$ be the probability distribution that a molecule in Brownian motion is at $\mathbf{x}$ at time $t$ and has not yet exited a domain $\omega$. If $\mathbf{x}_{0}$ is the starting position of the molecule diffusing with $\gamma_{0}$, then $c(\mathbf{x},t)$ fulfills

	$\displaystyle c_{t}(\mathbf{x},t)$	$\displaystyle=\gamma_{0}\Delta c(\mathbf{x},t),$	$\displaystyle\mathbf{x}$	$\displaystyle\in\omega,$		(3.1)
	$\displaystyle c(\mathbf{x},t)$	$\displaystyle=0,$	$\displaystyle\mathbf{x}$	$\displaystyle\in\partial\omega,$
	$\displaystyle c(\mathbf{x},0)$	$\displaystyle=\delta_{\mathbf{x}_{0}}.$

The homogeneous Dirichlet boundary condition here models that the particle is removed once it hits the boundary. The survival probability of the particle inside $\omega$ until time $t$ is then

$\displaystyle S(t)$

$\displaystyle=$

$\displaystyle\int_{\omega}c(\mathbf{x},t)d\omega.$

(3.2)

By Gauss’ formula the probability density $p_{\omega}(t)$ that the particle leaves $\omega$ at $t$ is

$$p_{\omega}(t)=-\frac{\partial S(t)}{\partial t}=-\gamma_{0}\int_{\partial\omega}\mathbf{n}\cdot\nabla c(\mathbf{x},t)ds,$$

(3.3)

where $\mathbf{n}$ is the outward normal. The expected time $E$ for the molecule to leave $\omega$ for the first time is given by

$$E=\int_{0}^{\infty}tp_{\omega}(t)dt=\int_{0}^{\infty}S(t)dt.$$

(3.4)

We now use this FET approach to compute the jump propensities $\lambda_{ii}$ in the space discrete mesoscopic model. We use a Cartesian grid with space discretization $h$ and $N$ nodes $\mathbf{x}_{i}$ in the domain $\Omega$. The voxels $\mathcal{V}_{i}$ are here defined by the dual mesh, see Fig. 3.1(b). Since the molecules are considered well mixed inside the voxels, the domain $\omega$ that diffusing particles have to leave to be well mixed in the next voxel has to include the centers of the neighboring voxels. On a Cartesian grid we showed in [47] that solving (3.1) on the circle $\omega_{i}$ with center $\mathbf{x}_{i}$ and radius $h$ (similarly a line of length $2h$ in 1D or a sphere with radius $h$ in 3D) and choosing $\mathbf{x}_{0}=\mathbf{x}_{i}$ gives the correct exit time from node $i$, see Fig. 3.1(b). Observe that $\omega_{i}\supsetneq\mathcal{V}_{i}$ and $\omega_{i}\cap\omega_{j}\neq\emptyset$ for neighboring nodes $i$ and $j$. Using (3.2) and (3.4) the expected exit time $E_{i}$ from $\mathcal{V}_{i}$ is

$\displaystyle E_{i}$

$\displaystyle=$

$\displaystyle\frac{h^{2}}{2d\gamma_{0}},$

(3.5)

see [23]. Since the jump propensity is the inverse of the exit time this agrees with the mesoscopic rate on Cartesian grids

$$\lambda_{ii}=\frac{2d\gamma_{0}}{h^{2}}=E_{i}^{-1}.$$

(3.6)

The probability to leave $\omega_{i}$ through a given part of the boundary $\partial\omega_{ij}$ at time $t$ is given by the proportion of fluxes

$\displaystyle\theta_{ij}(t)$

$\displaystyle=$

$\displaystyle\frac{\int_{\partial\omega_{ij}}\mathbf{n}\cdot\nabla c(\mathbf{x},t)ds}{\int_{\partial\omega_{i}}\mathbf{n}\cdot\nabla c(\mathbf{x},t)ds},$

(3.7)

and we can compute the expected probability to jump to a certain neighboring voxel by

$$\theta_{ij}=\int_{0}^{\infty}\theta_{ij}(t)p_{\omega_{i}}(t)dt=-\gamma_{0}\int_{0}^{t}\int_{\partial\omega_{ij}}\mathbf{n}\cdot\nabla c(\mathbf{x},t)dsdt.$$

(3.8)

Choosing $\partial\omega_{ij}$ to be the quarter segment of the boundary closest to $\mathbf{x}_{j}$, see the blue line in Fig. 3.1(b), yields the splitting probability $\theta_{ij}=0.25$ as expected on a Cartesian grid. The method has been extended to a rectangular grid and possible jumps to the diagonal neighbors in [55].

By conditioning on the first step, these time dependent equations can be converted to equations describing directly the expected quantities [65]. The expected exit time for a molecule starting to diffuse in $\mathbf{x}$ from the domain $\omega_{i}$ fulfills the Poisson’s equation

	$\displaystyle\gamma_{0}\Delta E(\mathbf{x})$	$\displaystyle=-1,$	$\displaystyle\quad\mathbf{x}\in\omega_{i},$		(3.9)
	$\displaystyle E(\mathbf{x})$	$\displaystyle=0,$	$\displaystyle\quad\mathbf{x}\in\partial\omega_{i}.$

Equivalently, the expected splitting probability for this molecule to exit through the boundary $\partial\omega_{ij}$ can be computed by the harmonic measure [60, Ch. 7], and fulfills the Laplace equation

	$\displaystyle\Delta\theta_{ij}(\mathbf{x})$	$\displaystyle=0,$	$\displaystyle\quad\mathbf{x}\in\omega,$		(3.10)
	$\displaystyle\theta_{ij}(\mathbf{x})$	$\displaystyle=1,$	$\displaystyle\quad\mathbf{x}\in\partial\omega_{ij},$
	$\displaystyle\theta_{ij}(\mathbf{x})$	$\displaystyle=0,$	$\displaystyle\quad\mathbf{x}\in\partial\omega_{i}\setminus\partial\omega_{ij}.$

We solve equations (3.9) and (3.10) and evaluate them at $\mathbf{x}_{i}$ instead of solving the time-dependent equations and the integrals above. It is sufficient to solve (3.10) three times for each node since $\sum_{j=1}^{4}\theta_{ij}=1$.

In the following, we will not simply use the circle $\omega_{i}$ to compute the mesoscopic rates, but we will prohibit the molecule from diffusing where the crowders are located.

The crowding molecules are represented as obstacles or holes in the domain $\omega_{i}$ with reflecting boundary conditions. Equations (3.1),(3.9) and (3.10) describe the diffusion of point particles. To account for the volume of the diffusing molecule its radius is added to the excluded volume for the center of mass, see Fig 3.1(a). We depict the crowders as circles with radius $R$, but it is important to mention that any shape is possible for the crowding molecules. The shape of the small (as compared to crowders) diffusing molecule is, however, restricted to circles or spheres with radius $r$.

Solving (3.9) and (3.10) numerically on the perforated domain $\omega_{i}$ means that the crowding molecules have to be resolved by a fine mesh. But, we have divided the global problem into $N$ local subproblems (one protective domain $\omega_{i}$ around each node $i$), which is only solved once and can be parallelized. This is similar to the approach in [7] where deterministic local equations are solved on media with porous microstructures. The stochastic simulation of the spatial SSA is then performed on the coarse mesh with $N$ nodes no longer resolving the individual obstacles. The boundary conditions on the global domain $\Omega$ (reflective or absorbing) are implemented by posing these conditions on the secants of the half or quarter circles, which are the protective domains for boundary nodes, see Fig. 3.2(a). In Fig. 3.2(b) we briefly illustrate how the first exit time approach can be further used to compute the jump rates for a Cartesian grid, discretizing a domain $\Omega$ with a curved boundary. Here the part of $\partial\omega_{i}$ originating from the circle is imposed with the boundary conditions in (3.9) and (3.10) and the part originating from $\partial\Omega$ with the boundary condition valid on $\Omega$, which usually are reflecting or (partially) absorbing. The circle is then divided in the same way as before to compute the splitting probabilities to the neighboring nodes remaining inside $\Omega$.

The simultaneous interpretation of the moving molecules being well-mixed inside the voxels $\mathcal{V}_{i}$ and jumping from node $\mathbf{x}_{i}$ to $\mathbf{x}_{j}$ leads to problems when including crowders. Consider the case where just the center $\mathbf{x}_{i}$ is blocked, but voxel $\mathcal{V}_{i}$ is sufficiently empty to be traversed, see Fig. 3.2(c). In this case the jump into voxel $\mathcal{V}_{i}$, is possible ($\lambda_{ji}>0$), but the expected time to leave $\mathcal{V}_{i}$ is infinity and hence the molecules get trapped inside $\mathcal{V}_{i}$. This does not agree with the microscopic situation, where molecules diffuse around $\mathbf{x}_{i}$. To avoid unrealistic trapping, all jump propensities $\lambda_{ji}$ to voxels whose vertex is isolated are set to zero. Equations (3.9) and (3.10) cannot be evaluated in $\mathbf{x}_{i}$ if the node is covered by the excluded volume and setting all $\lambda_{ji}$ to zero would over estimate the crowding effect considerably so we distribute the crowders such that $\mathbf{x}_{i}$ remains inside $\omega_{i}$.

By using the expected exit time from $\omega_{i}$ around $\mathbf{x}_{i}$ to calculate the jump coefficients in a crowded environment it is the crowder distribution inside the whole circle $\omega_{i}$ that affects the coefficients $\lambda_{ii}$ and $\theta_{ij}$. This differs from other approaches to simulate diffusion in a crowded environment with a discrete space jump process. In [45, 67, 75] it is only the percentage of occupied volume in the target voxel $\mathcal{V}_{j}$ and in [31] the difference in occupancy between $\mathcal{V}_{j}$ and the origin $\mathcal{V}_{i}$ that affect the jump rate $\lambda_{ij}=\theta_{ij}\lambda_{ii}$. In our approach the microscopic positions of all crowding molecules inside $\omega_{i}$ are resolved and influence the jump coefficients. In the case of non-spherical crowding molecules also the orientation is taken into account and long thin molecules with small volume can have a significant effect on $\lambda_{ij}$ and $\theta_{ij}$, see Fig. 3.3. Note that in contrast to normal diffusion the jump propensities are no longer symmetric, i.e. in general $\lambda_{ij}\neq\lambda_{ji}$ and $\lambda_{ij}\neq\lambda_{im}$ for $j\neq m$.

Solving $N$ local problems where complicated geometries have to be resolved, see Fig. 3.1(c) is computationally expensive and will be inefficient if the crowding molecules move and the coefficients $\lambda_{ii}$ and $\theta_{ij}$ have to be recomputed often. Since the crowders’ exact location is generally unknown we can compute statistics on a reference domain $\omega_{i}$ for given a percentage of occupied volume, a given shape and size of the molecules and $h$. Instead of solving $N$ PDEs of type (3.9) and $3N$ of type (3.10) at each time step, we can then sample the coefficients $\lambda_{ii}$ and $\theta_{ij}$ from these precomputed distributions. This will be especially applicable for moving crowders, where new coefficients can be drawn from the distributions on the time scale of their diffusion.

We first investigate how the jump propensity $\lambda_{ii}$ changes in different crowding situations. We therefore compute $\lambda_{ii}$ on a reference domain $\omega_{i}$ with $h=1$ and different crowders and sizes of the moving molecule $r$ and compare it to the jump rate $\lambda_{0,ij}$ in dilute medium. In Fig. 6.1, we compare the mean value of the jump propensities $\mathbb{E}[\lambda_{ii}]$ for different distributions of crowders and an increasing percentage of occupied volume $\phi$ with the jump propensities when no crowders are present, where $\mathbb{E}[\lambda_{0,ii}]=\lambda_{0,ii}=4$. We test two different sizes of crowding molecules for both rectangles and spheres. The reference line is the linear scaling where $\lambda_{ii}=(1-\phi)\lambda_{0,ii}$ as in [45, 63].

We observe that small obstacles (blue and orange in Fig. 6.1) hinder diffusion more than big crowders for the same percentage of occupied volume, since they have more reflecting surfaces than larger crowders. The same holds for elongated rectangular crowders (green and orange in Fig. 6.1), as they create long barriers without occupying a lot of volume. An increasing size $r$ of the diffusing molecule leads to an increase in the crowding effect, which is intuitive, since a bigger molecule finds less holes through which to escape. These results agree with the findings in [59] and [16]. We see that the averaged linear reduction of the jump propensity can be a good model when the diffusing molecule is about a tenth of the size of the crowders, but over- or underestimates the effect of occupied volume when the diffusing species is smaller or bigger, respectively. Since an average protein has a radius of ca. 2nm [63] and the biggest macromolecules in the cell, the ribosomes, have a radius of up to 15nm the linear correction is a good approximation for many scenarios. The reduction in jump propensity, however, starts to behave exponentially, as in [31], for large diffusing molecules. The case $r=0$ corresponds to a point particle which is irrelevant when simulating excluded volume effects, but we include it to show the limit for very small particles. To confirm these mesoscopic jump rates we compare them to the inverse of the expected exit time computed by a Brownian dynamics simulation. We simulate $1e4$ trajectories with the open source software Smoldyn [1, 2] for 10 different crowder distributions and equally sized crowding and moving molecules with $R=r=1e-1$ and see in Fig. 6.1(e) that the computationally expensive microscopic results agree well with the mesoscopic coefficients.

In Fig. 6.2, the standard deviation $\sigma(\lambda_{ii}/\lambda_{0,ii})$ initially increases as more and more crowders are added but converges towards zero when the system approaches the state where no escape to the boundary is possible.

Simply rescaling the jump propensity $\lambda_{ii}$ for each node $i$ by a constant factor will lead to normal diffusion at reduced rate. To observe anomalous diffusion in the crowded environment it is helpful to investigate the mean square displacement (MSD).

As mentioned in Section 2.2 the MSD is linear in time for normal diffusion

$$\langle\mathbf{x}^{2}(t)\rangle=2d\gamma_{0}t,$$

(6.1)

but for anomalous diffusion the relation is no longer linear

$$\langle\mathbf{x}^{2}(t)\rangle=2d\gamma_{0}t^{\alpha},$$

(6.2)

where $\alpha<1$ for subdiffusion. In [11, 57] it was shown that diffusion in a crowded environment can be modeled by a temporal change of the diffusion constant. First the molecules diffuses normally with rate $\gamma_{0}$ for very short time-scales, before it undergoes a transient anomalous phase with a changing diffusion coefficient $\gamma$ and finally stagnates into normal diffusion at a lower diffusion rate $\gamma_{\infty}$. The initial normal diffusion represents the time the molecule diffuses in the solution before it encounters the first adjacent macromolecule and is slowed down by collisions. On a large time scale the molecule appears to diffuse in a denser medium instead of around obstacles, hence the reduced diffusion rate $\gamma_{\infty}$, see Fig. 6.3. If the crowding macromolecules are distributed evenly the MSD decays monotonically between $\gamma_{0}$ and $\gamma_{\infty}$ (pale line), but due to stochastic variations in the medium it fluctuates before converging to $\gamma_{\infty}$ (dark line). In Fig. 6.3(a) we only depict collisions with the macromolecules responsible for excluded volume effects, but note that there are many collisions with the much smaller solvent molecules responsible for the Brownian motion.

Since we choose $h>R$ a jump in the mesoscopic model spans a number of macromolecules and the initial free diffusion phase is not resolved and we start to observe diffusion after the first jump of length $h$, that occurs after a critical time $t_{c}$ which can be approximated by

$$t_{c}\sim\frac{h^{2}}{2d\gamma_{0}}.$$

(6.3)

In the following we will plot the $\langle\mathbf{x}^{2}\rangle/2d\gamma_{0}t$ in log-log-scale for different crowding situations to examine when anomalous behavior occurs. Let $\mathbf{p}(t)\in\mathbb{R}^{N}$ be the probability vector for a diffusing molecule, such that $p_{i}(t)$ is the probability that the molecule is in voxel $\mathcal{V}_{i}$ at time $t$. As described in Section 2, $\mathbf{p}(t)$ evolves by the master equation

$$\mathbf{p}_{t}=\mathbf{D}\mathbf{p}(t),\quad\mathbf{p}(0)=\mathbf{p}_{0},$$

(6.4)

where $D_{ij}=\lambda_{ji}$ for $i\neq j$ and $D_{ii}=-\lambda_{ii}$. The initial probability distribution $\mathbf{p}_{0}$ is $(0,\dots,1,\dots,0)^{T}$ with one at the starting node $\mathbf{x}_{0}$. We choose the discretization such that $\mathbf{D}$ is small enough to solve (6.4) numerically and compute the mean square displacement by

$$\langle(\mathbf{x}(t)-\mathbf{x}_{0})^{2}\rangle=\sum_{i=1}^{N}p_{i}(t)(\mathbf{x}_{i}-\mathbf{x}_{0})^{2}.$$

(6.5)

In the following experiments we discretize the square $[0,1]\times[0,1]$ into $N=n^{2}$ voxels with space discretization $h=1/(n-1)$. If not mentioned otherwise we release the molecules in $[0.5,0.5]$ at time $t=0$ and choose $n=41$. To avoid boundary effects we set homogeneous Dirichlet boundary conditions on $\partial\Omega$ and show the solutions as long as more than $99\%$ of the mass is preserved, i.e. $\sum_{i=1}^{N}p_{i}(t)>0.99$.

The MSD is only one quantity of interest to examine, but since it is a mean not all features are captured and we will now compare the distributions of molecules resulting from either a mesoscopic or macroscopic simulation. Again, we discretize the square $[0,1]\times[0,1]$ with homogeneous Dirichlet boundary conditions into 41 nodes in each direction and let molecules start diffusing in $(0.5,0.5)$ in an environment with rectangles of different sizes. We solve (6.4) once with the mesoscopic $\mathbf{D}$ and once with $\tilde{\mathbf{D}}$, where the off-diagonal elements are all equal to $\lambda_{ii}/4$, corresponding to a finite difference approximation of the macroscopic equation with the space dependent diffusion constant $\gamma(\mathbf{x})$ derived from the $\lambda_{ii}$s. In Fig. 6.5 the macroscopic model agrees with the mesoscopic results for small and evenly distributed crowding molecules , whereas for long barriers only the mesoscopic approach simulates the expected behavior. This is due to the symmetrization of $\tilde{\mathbf{D}}$, so that only $\mathbf{D}$ can capture the asymmetric diffusion close to the barriers. The diagonal barriers are not completely impermeable in the mesoscopic model since a small part of the boundary $\partial\omega_{ij}$ remains.

Due to the reduction of diffusion in a crowded environment and (5.2) the overall reaction rate is decreased. It has, however, been shown [28, 68], that protein associations can also be enhanced. We examine the mean time $\mathbb{E}^{r}_{j}$ until the bimolecular reaction (6.6)

$$A+B\xrightarrow{k_{r}}C$$

(6.6)

happens, where reactant $A$ is confined to voxel $\mathcal{V}_{i}$, and molecule $B$ starts diffusing in voxel $\mathcal{V}_{j}$ at time $t=0$, see Fig. 6.6(a). Due to molecular crowding we assume a simplified space dependent diffusion map with $\gamma_{i}<\gamma_{0}$ inside $\mathcal{V}_{i}$ and $\gamma<\gamma_{0}$ in the rest of the domain. With $k_{i}$ given by (5.2) and $\lambda_{ii}$ by (3.6) we can use conditioning on the first step to compute the expected time until the reaction happens:

	$\displaystyle\mathbb{E}^{r}_{i}$	$\displaystyle=$	$\displaystyle\frac{k_{i}}{k_{i}+\lambda_{ii}}\frac{1}{k_{i}+\lambda_{ii}}+\frac{\lambda_{ii}}{k_{i}+\lambda_{ii}}\left[\frac{1}{k_{i}+\lambda_{ii}}+\sum_{m}\theta_{im}E(\mathbf{x}_{m})+\mathbb{E}^{r}_{i}\right]$		(6.7)
	$\displaystyle\mathbb{E}^{r}_{j}$	$\displaystyle=$	$\displaystyle E(\mathbf{x}_{j})+\frac{h^{2}+4\gamma_{i}\sum_{m}\theta_{im}E(\mathbf{x}_{m})}{h^{2}k_{i}},$		(6.8)

where $E(\mathbf{x}_{m})$ is the expected time it takes a molecule located in the neighboring voxel $\mathcal{V}_{m}$ to jump into voxel $\mathcal{V}_{i}$ and can be computed by solving

	$\displaystyle\mathbf{p}_{t}(t)$	$\displaystyle=$	$\displaystyle(\mathbf{D}-\mathbf{K}_{i})\mathbf{p}(t),\quad\mathbf{p}(0)=\mathbf{p}_{0},$		(6.9)
	$\displaystyle E(\mathbf{x}_{m})$	$\displaystyle=$	$\displaystyle\int_{0}^{\infty}\sum_{k=1}^{N}|\mathcal{V}_{k}|p_{k}(t)dt,$		(6.10)

where $\mathbf{K}_{i}$ models a sink at node $\mathbf{x}_{i}$ and is the zero matrix except for $K_{i,i}=10^{9}$ and $\mathbf{p}_{0}$ is the zero vector except for $p_{0,m}=1/h^{3}$, see [54] for a derivation. We solve these equations numerically in 3D for a cube with length $L=1$ and a uniform discretization with $h=0.1$ in space and reflecting boundary conditions. The voxel $\mathcal{V}_{i}$, where the reaction happens is chosen to be the center voxel, such that $E(\mathbf{x}_{m})$ are equal for all 6 neighbors. The diffusing molecule $B$ starts in $\mathbf{x}_{j}=(0.7,0.5,0.5)$. In Fig. 6.6 we compare the mean binding time in the crowded $\mathbb{E}_{j}^{r}$ environment with different $\gamma_{i}$ and $\gamma$ to the time $\mathbb{E}^{r}_{j,0}$ it takes to react in a dilute solution where $\gamma(\mathbf{x})=\gamma_{0}=1$. The data points with scaled error bars ($\pm(\frac{\sigma}{\sqrt{M}\mathbb{E}^{r}_{j}}+\frac{\sigma_{0}}{\sqrt{M}\mathbb{E}^{r}_{j,0}})$) are from a SSA simulation of the reaction-diffusion process with $M=200$ trajectories for $k_{r}=1e-4$ and $k_{r}=1e-3$ and $M=500$ for $k_{r}=5e-3$ and $k_{r}=1e-2$.

An overall slower diffusion rate $\gamma<\gamma_{0}$ as a result of obstacles reduces the rate of bimoleculear reactions (5.2) in each voxel. But, due to an uneven distribution of crowding agents compartmentalization with locally differing diffusion rates can occur inside the cells. In Fig. 6.6 the overall reaction rate can be increased for a compartmentalization where $\gamma_{i}$ is smaller than $\gamma$, despite the locally slower reaction. Because once the diffusing molecule enters $\mathcal{V}_{i}$ it escapes slower and hence gets trapped close to its reaction partner, which increases the chance of collisions. Like this, cells can boost their efficiency by locating important reaction complexes in areas of slower diffusion, which will be especially productive for reaction cascades where intermediate products are already produced inside the compartment and have a low chance of escaping before being processed further. In the limit when $\gamma_{i}\to 0$ the binding time goes towards infinity, since then both reactants are immobile inside $\mathcal{V}_{i}$ and never collide.