Diffusion through thin membranes : Modeling across scales

All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Diffusion through thin membranes: Modeling across scales Aho, Vesa; Mattila, Keijo; Kühn, Thomas; Kekäläinen, Pekka; Pulkkinen, Otto; Minussi, Roberta Brondani; Vihinen-Ranta, Maija; Timonen, Jussi

From macroscopic to microscopic scales it is demonstrated that diffusion through membranes can be modeled using specific boundary conditions across them. The membranes are here considered thin in comparison to the overall size of the system. In a macroscopic scale the membrane is introduced as a transmission boundary condition, which enables an effective modeling of systems that involve multiple scales. In a mesoscopic scale, a numerical lattice-Boltzmann scheme with a partial-bounceback condition at the membrane is proposed and analyzed. It is shown that this mesoscopic approach provides a consistent approximation of the transmission boundary condition. Furthermore, analysis of the mesoscopic scheme gives rise to an expression for the permeability of a thin membrane as a function of a mesoscopic transmission parameter. In a microscopic model, the mean waiting time for a passage of a particle through the membrane is in accordance with this permeability. Numerical results computed with the mesoscopic scheme are then compared successfully with analytical solutions derived in a macroscopic scale, and the membrane model introduced here is used to simulate diffusive transport between the cell nucleus and cytoplasm through the nuclear envelope in a realistic cell model based on fluorescence microscopy data. By comparing the simulated fluorophore transport to the experimental one, we determine the permeability of the nuclear envelope of HeLa cells to enhanced yellow fluorescent protein.

I. INTRODUCTION
Membranes are selective barriers restricting transport of material between two regions. They are ubiquitous in living organisms; the transport of molecules between a cell and its exterior, and between many cellular compartments, is regulated by lipid bilayers outlining them. These membranes are extremely thin in comparison with the size of the whole cell. For example, the nuclear envelope of a cell is about tens of nanometers thick [1], while a typical eukaryotic cell is of the membrane. In other words, the membrane induces a concentration jump δρ across it. Moreover, the coefficient of proportionality P is the membrane permeability. In the standard solubility-diffusion model, for example, P = KD m /w, where K is a dimensionless partition coefficient, D m is the intrinsic diffusion coefficient in the membrane, and w is the membrane thickness [4,5].
In this article, we examine diffusion through thin membranes. In a macroscopic analysis, we treat the membrane as a transmission boundary condition between two regions, and derive solutions for the concentration fields in simple configurations. However, many applications of great interest are beyond an analytical treatment due to inherent complexities involved. Hence it is of practical importance to consider numerical modeling of membranes as well. A numerical treatment of continuum-based models requires discretization of the space. When the membrane is considered as a singularity also in the scale of discretization, say in the scale of lattice spacing, it is not immediately clear how the permeability in a numerical scheme should be implemented to be in accordance with the permeability of Eq. (1) in the continuum description. The situation is even more complex in numerical methods which do not explicitly contain the macroscopic differential equations they solve, such as the lattice-Boltzmann method (LBM) [6,7].
In order to focus on this issue, we discretize the macroscopic transmission conditions so as to derive an expression for the mass flux between two discretization points separated by a membrane. The obtained expressions can be readily utilized, e.g., in finite-difference and finite-volume methods for solving the diffusion equation. In a mesoscopic scale, we propose an implementation of LBM, where the membrane is treated as a partial-bounceback condition between two lattice nodes. After presenting the discrete scheme, we perform a multiple-scale analysis and derive an explicit expression for the mass flux at the membrane. To bridge the gap between microscopic and mesoscopic descriptions, we show, in a microscopic description, that the flux coincides with the inverse mean waiting time of a single molecule in the vicinity of the membrane.
Finally, we show that the numerical approximations with LBM agree with analytically derived solutions at a macroscopic level, and show the applicability of the LBM scheme by determining the permeability of the nuclear envelope to enhanced yellow fluorescent protein (EYFP). In this numerical determination of permeability, photobleaching experiments are compared against simulations executed using realistic, three-dimensional, image-based geometries.

II. CONTINUUM THEORY
To begin with, in a macroscopic scale the diffusion equation with a constant diffusion coefficient D is given by and Fick's law for the mass flux J is given by Using Eq. (1) and continuity of the mass flux, we find that the transmission conditions at the membrane (at x = 0) can be , illustrates an interpretation of the effective membrane thickness, 2 , using = D/2P . Shifting the concentration profile by will remove the concentration jump and, hence, remedy the discontinuity [8]. The bottom part of the figure, (c) and (d), is a schematic description for the linear interpolation of ρ(0 + ) and ρ(0 − ) according to Eq. (11). The stars indicate the two interpolation points. expressed in the form, where ρ(0 + ) and ρ(0 − ) now denote concentrations at the opposite faces of the membrane. Note that the mass flux across the membrane, Eq. (1), can also be expressed in the form, an effective thickness of the membrane, 2 , is defined using = D/2P [8]. That is, even though the concentration is not formally differentiable across the membrane, a diffusive transport across it can be expressed by the usual Fick's law, when the partial derivative is understood in the manner described above (see the top part of Fig. 1).

A. A steady-state configuration
The steady-state solution for a one-dimensional (1D) diffusion equation, given by J = −D∂ x ρ = constant, implies a linear concentration profile. Let us consider a domain [−L/2,L/2] with a given reference concentration ρ r , and a fixed concentration difference ρ L between the domain ends, i.e., the boundary conditions, can support a nonzero, constant mass flux at steady state. Since the configuration is symmetric, A transmission boundary condition then enforces a continuous mass flux for a discontinuous concentration, and can be 043309-2 expressed in the form, Since the concentration profile in the bulk is linear at steady state, Exactly the same expression is obtained for (∂ x ρ)(0 + ). Hence, which finally leads to where γ = LP /D.

B. A transient configuration
An analytical solution can be derived for the above membrane boundary condition in the case, where x ∈ [−L,0) ∪ (0,L], and concentration is initially uniform on one side of the membrane and zero on the other. Mathematically, we write a no-flux boundary condition at ±L, and an initial condition at t = 0: Note that in the following analysis we will, for simplicity, use dimensionless densities which can be understood as normalized quantities of the formρ = ρ/ρ r , where ρ r is a reference density. Defining ξ = x/L andτ = t D/L 2 , we can solve the resulting dimensionless problem by separation of variables (see Appendix A), and find that where again γ = LP /D, and the eigenvalues μ satisfy

III. DISCRETIZATION OF THE TRANSMISSION BOUNDARY CONDITIONS
In a fine enough discrete numerical model, the membrane can be constructed explicitly by assigning a number of discretization points to it, and setting the membrane diffusion coefficient D m in those points such that the desired flux through the membrane, Eq. (3), will be obtained. Here we continue instead to treat the membrane as a boundary condition. In a true multiscale-modeling fashion, we place the membrane between two discretization points and discretize the transmission conditions of Eq. (4). We set up our computational grid such that there is a grid point on each side of the membrane with an equal distance h to it. Furthermore, the continuity condition for the flux at the membrane, approximated to first order in space with backward and forward differences, is given by Setting P [ρ(0 + ) − ρ(0 − )] equal to both sides of Eq. (10), separately, we can express ρ(0 + ) and ρ(0 − ) in terms of ρ(h) and ρ(−h). In this way we find that That is, Eq. (11) defines a linear interpolation for obtaining the concentrations at the membrane using the two adjacent discretization points (see the bottom part of Fig. 1). Note that this interpolation is consistent, when h → 0. After substituting these expressions into Eq. (1), we find further that where ρ = ρ(h) − ρ(−h) and the dimensionless coefficient λ = h/ = 2hP /D. That is, the mass flux across the membrane is now expressed using the concentration difference between the two adjacent grid points, not between the two faces of the membrane. Therefore the measured difference must be adjusted so as to be consistent with the definition of the permeability coefficient. Thus ρ/(1 + λ) can be considered as the concentration jump interpolated at the membrane.

A. Effective diffusion coefficient in the membrane
An alternative interpretation of the mass flux in Eq. (12) can be adopted. Using the definition of λ, we arrive at This same expression, but in a different context, has been reported, e.g., in Refs. [9,10]. Now the partial derivative at the membrane is approximated using the discretization points and, in particular, not across the effective membrane thickness, cf. Eq. (5). While the coefficients in Eq. (13) originate from the concentration-jump approximation, this mass-flux expression is quite elegant: The effect of the membrane is manifested by a single coefficient λ, having the role of reducing the corresponding bulk mass flux. Accordingly, the coefficient λD/(1 + λ) can be considered as the effective diffusion coefficient for the membrane. The obtained expressions, Eqs. (11)- (13), can be utilized, e.g., in finite-difference and finite-volume methods for solving the diffusion equation. 043309-3

IV. LATTICE-BOLTZMANN MODELING OF DIFFUSION ACROSS A MEMBRANE
Numerical modeling on the level of distribution functions of the solute, i.e., on a mesoscopic scale, is conveniently implemented using the lattice-Boltzmann method [6,7], which has emerged as an alternative in a computational treatment of complex transport phenomena. Diffusion processes, for example, can be simulated using the standard lattice-Boltzmann evolution equation with the BGK collision operator [11], where q is the number of discrete velocities in the model, τ * is a dimensionless relaxation time, ρ denotes the local concentration at a given instant, and w i are suitable weight coefficients for the equilibrium function, f eq i = w i ρ. The lattice spacing δr, and discrete time step δt, define a reference velocity, c r = δr/δt, and c i = c r c * i , where c * i are the dimensionless lattice vectors. The first few moments of the discrete equilibrium distribution function are given by where b is a constant fixed by selection of the weight coefficients w i (see Appendix B). The diffusion coefficient can be derived, e.g., with the Chapman-Enskog procedure [11]: As already mentioned, the most elementary way to construct the membrane in a lattice-Boltzmann framework is to assign a number of lattice nodes to it; the flux is controlled with the diffusion coefficient D m , e.g., by setting the relaxation time appropriately. This is also the approach taken in Ref. [12] to model the nuclear envelope. However, if the membrane restricts the transport strongly, relaxation time at the membrane sites can be inconveniently small (with respect to the numerical stability [13]), or too large in the surrounding environment (with respect to the numerical error [11]). A possible way to circumvent this problem is to use a local grid refinement at the membrane and its immediate neighborhood. Unfortunately in a complex geometry, which is the typical case, such refinement may be cumbersome or not feasible at all.
Instead, we follow Ref. [14], where a mesoscopic boundary condition corresponding to Eq. (1) was reported, and utilize this simple boundary condition for the implementation of a membrane in a lattice-Boltzmann scheme. The condition relies on a single parameter (we will call it φ in the following) that controls the proportion of particles which are able to pass through the membrane. Namely, we introduce a partial halfway-bounceback rule for the distributions that cross the membrane, where subscript −i denotes the direction which is opposite to that of i, i.e., c −i = −c i . That is, the location of the membrane is halfway between two neighboring lattice nodes. Depending on the application, φ can be associated with, e.g., the relative area of the membrane, which allows diffusive transport. Equation (15) provides a subgrid-scale model for the membrane diffusion, and it is conveniently implemented even in arbitrarily complex geometries. Previously partialbounceback schemes have been used in the LBM context for modeling transport in porous media [15][16][17][18], thermal conduction across the boundary between two materials [19], or slip flows [20,21].
The Chapman-Enskog analysis of the scheme of Eq. (15) shows (see Appendix B) that it leads to the modified flux as obtained above at the membrane, Eq. (12), and the restrictive effect of the membrane is included in the coefficient λ = P δr/D, h = δr/2, where the permeability is given by That is, the partial-bounceback scheme is an alternative discretization of the transmission boundary condition. It is also evident that this expression for the permeability implies a reasonable limiting behavior, Furthermore, when D approaches infinity, i.e., when transport is limited only by the effect of the membrane, J m approaches the value −P δr ∂ α ρ. This is consistent with Eq.
(1), and shows that P can indeed be thought of as the permeability of the membrane. This allows us to compare an approximate solution computed by the lattice-Boltzmann scheme with the analytical solution derived above.

V. MICROSCOPIC MODELING OF DIFFUSION THROUGH A MEMBRANE
We interpret the permeability of the Chapman-Enskog analysis, Eq. (16), in terms of the kinetics of individual solute molecules by studying a single-particle stochastic diffusion process, the probability density of which coincides with the solution of the macroscopic diffusion equation, Eq. (2), with the flux boundary condition, Eq. (4), at the membrane. In that sense, it is the exact microscopic counterpart of the macroscopic description.
The idea of the microscopic process is that a particle must on average collide with the membrane a sufficient number of times to find a pore, through which it passes to the other side of the membrane. Hence there is a pore-density-dependent probability to go through the membrane given that the particle is at its immediate vicinity. In a continuum model, the time spent at the membrane is measured by the local time [22], where X t is the position of the particle at time t and I is an indicator function. The local time is clearly a stochastic process that can only increase in time. The trajectories of the process L t typically consist of long periods with no change in its value, as the Brownian particle is on an excursion far 043309-4 away from the membrane, and bursts of rapid increase when the particle is in its close vicinity.
A continuum stochastic process X t that passes through the membrane when its local time at the membrane exceeds an exponential random variable is an extension of the elastic Brownian motion [20], called snapping-out Brownian motion in a recent preprint [23]. In the preprint, the process was also shown to be the limit of a diffusion through a membrane of a vanishing diffusion constant D(w) as the membrane width w → 0. The first passage through the membrane was further shown to be marked by the local time exceeding an exponential random variable with a characteristic time scale, τ c = lim w→0 w/D(w), which depends on the microscopic structure of the membrane. However, it can be written down for the residence times of random walks with an idealized, infinitely thin membrane (see Appendix C).
Taking the continuum limit results in a snapping-out Brownian motion, and provides the equivalence of singlemolecule and ensemble-averaged permeabilities, In other words, the inverse mean local time at the membrane before the particle passes through to its other side is proportional to the permeability in an ensemble experiment. This result is reminiscent of single-molecule reaction kinetics, in which case the mean waiting time to the next chemical reaction that involves an individual enzyme molecule equals the inverse of reaction velocity in a large ensemble of constituent molecules [24].

VI. VERIFICATION OF THE LB MEMBRANE SCHEME
Here we numerically verify that the above presented LB membrane scheme indeed is a consistent discretization of the transmission boundary conditions presented in Sec. II. To this end, we simulated diffusion across a membrane in a one-dimensional setting using the standard LB scheme for diffusion [11] together with the D1Q3 discrete velocity set. We chose b = 2/3 which implies w i = 1/3 for all i. Furthermore, we used a constant value for the dimensionless relaxation time, i.e., τ * = τ/δt = 3/2. The permeability was controlled with φ according to Eq. (16). The values of D and P used in the simulations were then substituted into Eqs. (7) and (8) when evaluating the analytical profile.

A. Steady-state benchmark case
First we considered a steady-state benchmark case. The main idea is to compare the analytically derived approximation for the mass flux, Eq. (12), with the mass flux computed numerically. We also compare the analytical and simulated concentration profiles.
In a standard LB scheme for diffusion, the local mass flux is computed as the first moment of nonequilibrium distributions averaged over pre-and post-collisional states. Let J n denote the flux across the membrane (where n stands for the numerical flux). We then measure (when using D1Q3) FIG. 2. The 1D steady-state benchmark configuration. The LB scheme is utilized together with the D1Q3 discrete velocity set. The two nodes adjacent to the membrane are located at −δr/2 and δr/2.
where f neq E (δr/2) is the nonequilibrium distribution at the right-hand side of the membrane and at the pre-collision state; f neq W (−δr/2) is defined accordingly (see Fig. 2).
We carried out the simulations enforcing constant concentrations ρ(−L/2) = ρ r + ρ L /2 and ρ(L/2) = ρ r − ρ L /2 at the boundaries: ρ r and L are the reference concentration and domain length, and ρ L is the given concentration difference across the domain. These boundary concentrations were also used as the initial concentrations on the left and right subdomains.
All tests were made using N = 40 lattice nodes. When comparing to an analytical solution with given P , D, and L, the simulation parameters are then the lattice spacing We performed two groups of simulations. With the first group we compared the steady-state concentration profiles obtained numerically with the analytical solution, Eqs. (6) and (7), for two different values of P and three different values of ρ L . With the second group of simulations we compared the mass fluxes obtained from Eqs. (12) and (17) using three different values for ρ L and ten values of P for each ρ L . Figures 3 and 4 show the results of these simulations. From the figures it is apparent that there is an excellent agreement between the numerical and analytical solutions. In fact, the analytical solutions were reproduced with machine precision, when comparing both the concentration profiles and the mass fluxes. These results are not too surprising as the presented LB scheme for modeling diffusion across membranes is secondorder accurate, based on theoretical considerations, and the benchmark is a linear configuration.

B. Transient benchmark case
Next we considered a transient benchmark case for verifying the LB membrane scheme. When comparing the analytical solution, Eq. (8), with the simulated solution, we represented the simulation domain with 100 lattice nodes and the membrane was again located at the center of the domain. The initial concentration was ρ 0 = ρ r on the left side of the membrane and ρ 0 = 0 on its right side. At the borders of the simulation domain the standard bounceback rule [25] was applied (i.e., a noflux condition). The values of D and L were fixed to the values used in the previous section, and the values of P and the time t were varied. The number of grid points determines the lattice spacing (L is fixed) as well as the discrete time step δt due to fixed D and τ * . The permeability parameter φ and the number of simulation time steps were then adjusted to deliver the given P and t. The results of these comparisons are presented in Fig. 5. In each case we observed a very good agreement between the simulated and analytically determined solutions.
To numerically confirm the second-order accuracy, predicted analytically by using the Chapman-Enskog analysis, we measured the convergence rate of error with respect to the grid spacing δr. We repeated the transient simulation and varied the number of grid points in the simulations from 10 to 500. The (absolute) error was measured by comparing the analytical and simulated density values at the lattice point nearest to the membrane on the left side of it. From this experiment, we observed that the convergence rate of error is of second order in space as expected (see Fig. 6).

VII. APPLICATION OF THE LB MEMBRANE SCHEME
After a successful verification of the LB membrane scheme, we demonstrate its capability as a computational tool in a realistic research problem. Namely, we apply the scheme in a biological context, and utilize it for determining the permeability of the nuclear envelope of HeLa cells to enhanced yellow fluorescent protein (EYFP).
Proteins can travel through the nuclear envelope passively by diffusion, or they can be actively transported by proteins called karyopherins. The karyopherin-mediated active transport requires that the transported protein contains peptide signals called nuclear localization signals or nuclear export signals in the case of nuclear import or nuclear export, respectively. The green fluorescent protein (GFP) and GFPderived fluorescent proteins, like EYFP used here, do not contain these signals and are known to be small enough to travel through the envelope passively.
The irregular shape of the cell and its nucleus, the nonuniformity of the equilibrium concentration of a molecule of interest, and movement of particles between cell organelles and their exterior are some of the problems that need to be addressed when considering the diffusion of particles within cells. These complexities usually mean that analytical methods are very difficult to utilize, and numerical methods become necessary. Instead of using a microscopic model for diffusion, like a random walker to model individual particles, it is more useful to work here on the particle density level instead, since in our approach we compare our simulations to experimentally determined densities. The connection between the particle densities and the lattice-Boltzmann distribution functions is simple (see Sec. IV) making it easy to interpret our simulations. Additionally, Eq. (16) gives us access to the permeability value from the simulation parameters utilized at the location of the nuclear envelope.
With these considerations in mind, we determine the permeability of the nuclear envelope to EYFP by comparing the measured evolution of the EYFP distribution in the cells to the evolution simulated with the lattice-Boltzmann method using different values of permeability. That is, we simulate nuclear transport in a real cellular geometry. The realistic, image-based simulation geometry is itself obtained by using fluorescence microscopy. By correlating the experimental and computational results, it is then possible to determine a permeability value for the nuclear envelope from the simulation parameters utilized.

Cells and culturing
Human-carcinoma HeLa MZ cells were used in all the experiments. The cells were grown in Dulbecco's modified Eagle medium with 10% FBS, while incubated at 37°C temperature in the presence of 5% CO 2 . The cells were transfected with plasmids encoding a freely diffusing yellow fluorescent protein (EYFP-N3) and a histone-bound cyan fluorescent protein (H2B-ECFP) using the TransIT-2020 reagent (Mirus Bio LLC, Madison, WI) in 5-cm glass-bottom dishes (MatTek, Ashland, MA) 24 hours before the measurements.

Microscopy
An Olympus FV1000 confocal microscope was used in the experiments with a 60× water-immersion objective whose numerical aperture was 1.2. The stage heater of the microscope was set to 37°C. The equilibrium distributions of EYFP and H2B-ECFP were imaged in the whole cell by taking confocal planes from the top to the bottom of the cell. The dynamics of EYFP was imaged in one confocal plane only. The voxel size was set to (150 nm) 3 in three-dimensional (3D) imaging and (150 nm) 2 in two-dimensional (2D) imaging. EYFP was excited using the 514-nm line of an argon laser, and H2B-ECFP using the 405-nm line of a diode laser.
The equilibrium distributions I eq of EYFP and H2B-ECFP were first imaged in the whole cell using an image size of 320 by 240 pixels. The average thickness of the cells was about 80 pixels. Since the EYFP distribution that is seen with the microscope is, macroscopically speaking, in equilibrium, we needed to disturb this equilibrium in order to measure its dynamics. This was achieved using the photobleaching technique, where a high-intensity laser pulse is directed to a region in a cell, destroying the fluorescence in that area. It was observed that after the photobleaching the diffusion of EYFP back towards the equilibrium was very fast, so the process could not be imaged in the entire cell but only in one confocal plane. To further increase the acquisition speed, a smaller area of 320 by 90-100 pixels was selected such that the area contained some of the nucleus, cytoplasm, and nuclear envelope. The dynamics of EYFP was then observed in this area by first photobleaching a circular area of 1.5-μm diameter in the nucleus using the 514-nm laser line with maximum power, and taking images with 5.2-5.5 frames per second thereafter. The imaged intensity after the bleaching in this one layer is called I exp in the following.

B. Simulation setup
For simulations the nucleus was segmented using fluorescence data of H2B-ECFP, and the whole cell was segmented from the background using the EYFP data [see Figs. 7(a)-7(d)]. As the initial condition for the simulation we selected the EYFP distribution after the photobleaching [Figs. 7(e) and 7(f)]. Even though we imaged the photobleaching process in one confocal plane, the laser beam causes the bleached profile to extend also in the direction perpendicular to the imaged plane. We assumed that the intensity profile of the laser was constant in the z direction [12] and calculated a fractional loss of intensity, γ (x,y), such that for the before-and after-bleach intensities the relation I after = γ I before held. The whole-cell initial concentration distribution for the simulation, I sim (r,t = 0), was then set by multiplying the intensities of the imaged equilibrium concentration of EYFP in every layer with γ (x,y) (extrapolation of the profile in the z direction).
The cell is filled with microscopic structures which prevent the fluorophores from entering certain regions of the cell, effectively reducing the accessible volume for the fluorophores. As these structures cannot be resolved by the microscope, the measured equilibrium distribution appears to be heterogeneous and nonuniform [see Fig. 7(a)]. Due to this nonuniformity, we had to solve the diffusion equation for the dimensionless (relative) intensity, I * (r,t = 0) = I sim (r,t = 0)/I eq (r) [26]. Solving the equation for I sim , instead of I * , would eventually lead to a uniform equilibrium concentration, which, as mentioned above, is not correct.
For these simulations we used the standard LB scheme for diffusion with the D3Q7 discrete velocity set; parameter b was fixed to 2/7. The grid spacing δr was chosen to be the same as the voxel size of the imaged cells (150 nm). At the borders 043309-7 of the cell a full bounceback rule [25] was applied (i.e., a no-flux condition), and the flux through the nuclear envelope was implemented by the membrane model presented in Sec. IV. This was achieved by applying the partial-bounceback rule of Eq. (15) to a distribution function f i (r) if r was located in the segmented nucleus (cytoplasm) and an adjacent point r + δrê i was located on the other side of the membrane in the segmented cytoplasm (nucleus). The permeability of the envelope was assumed to be uniform, since the density of the nuclear pores is high and they are rather uniformly located over the nuclear envelope [27,28]. Furthermore, the permeability was assumed to be constant during the simulation, since the total simulated time (only a few seconds) was short.
Besides the volume exclusion effect mentioned above, the cellular structures also reduce the mobility of molecules in their vicinity [12]. Therefore we had to scale the diffusion coefficients in the cell based on the amount of volume exclusion in the imaged equilibrium distribution. Values of τ (r) were set based on a simple capillary model, where ε is a porosity parameter that was assumed to be directly proportional to fluorescence intensity in the equilibrium fluorophore distribution, that is, ε(r) = I eq (r)/I max , where I max is the maximum fluorescence intensity of EYFP in the cell. For τ this means that we had to set them according to the equation, where τ * max is the maximum value of τ * in the cell, since, according to Eq. (14), D ∝ τ * − 1/2. The value of τ * max can now be chosen quite freely, e.g., based on accuracy, stability, or workload considerations for the problem. We found that τ * max = 2 was a good choice for our problem. Note that in our analytical considerations, we assumed that the diffusion coefficients on either side of the membrane are the same. This is strictly speaking not true in our simulations, but the gradient ∇D(r) is still sufficiently small to justify this assumption. By varying φ, the transmission parameter for the nuclear envelope, we obtained several diffusion patterns that were analyzed and fitted to experimentally measured diffusion to yield δt as a fitting parameter. In other words, φ was used  to control the amount of nuclear envelope permeation in the simulation, while δt was used to connect the whole simulation to the experiment to which it was compared. This works because both D and P are inversely proportional to δt. The fitting procedure is described in the next section.

Data analysis
Distributions of I * that were obtained from the simulations were compared to measured relative fluorescence of the corresponding experiments, I * exp = I exp /I eq , by calculating cross-correlation functions given by Cc(exp,sim) = where subscript "exp" refers to experiment and "sim" to simulation, I (x,y) are fluorescence intensities of the pixels, I aver average intensities of the images, σ s the standard deviations of the images, and N is the number of pixels in the image. Notice that the simulations were conducted in the whole cell, while the experimental EYFP diffusion was only imaged in one confocal plane of the cell. For this reason, the corresponding region had to be extracted from the simulations for the comparison. For each experimental time point, the simulation time point that gave the largest cross-correlation value between the two images was determined. In this way a curve was obtained, the linearity of which was used as an indicator of the accuracy of the simulation. For the simulation that yielded the most linear curve, the value of the simulation time step δt was extracted from the slope of the curve. The values of δt and φ were then used in Eq. (16) to determine the nuclear envelope permeability.

C. Determination of the nuclear permeability of EYFP
We determined the nuclear envelope permeability of EYFP in HeLa cells by comparing the evolution of EYFP distribution in the cells, measured by fluorescence microscopy, to evolution simulated with different values of permeability (see Fig. 8). Figure 9 shows that when the value of φ in the simulation was too high, the simulation progressed too fast compared to the experiment. Similarly, the simulation progressed too slowly when φ was too low. This allowed us to find an optimal FIG. 8. A comparison of the experimental (left-hand column) and simulated (right-hand column) intensities I * at three time points. Notice that the experimental data contains imaging noise. This is also true for the initial condition of the simulation, since it is determined from the experimental data. φ, and to determine the nuclear permeability as described above.
We performed simulations for 12 cells, and the determined permeabilities varied between 0.16 and 0.98 μm/s. The average permeability and its standard deviation were P = 0.5 ± 0.3 μm/s. The determined simulation time step, together with the used simulation parameters, allowed us to calculate also the diffusion coefficients in the cell. The determined average diffusion coefficients in the nucleus and cytoplasm were D aver nuc = 27 ± 7 μm 2 /s and D aver cyt = 15 ± 4 μm 2 /s. Using the solubility-diffusion model (see Introduction), P = KD m /w, together with equation A pore = K/n, where A pore is the cross-sectional area of one nuclear pore channel and n is the area density of the nuclear pores, we calculated a rough estimate for the effective nuclear pore diameter. The value D m that we used was the average of the nucleus and cytoplasm diffusion coefficients that we determined (21 μm 2 /s) and the value of n was the average of values from [27,28] (25 μm −2 ) and w ∼ 65 nm [29,30]. This way we obtained the effective nuclear pore area A pore ∼ 62 nm 2 , or effective pore diameter d pore ∼ 8.9 nm, which is in good agreement with experimental data [31,32].

D. Discussion
The average diffusion coefficients that we obtained in the nucleus and cytoplasm are close to what has been found in previous studies for EYFP and structurally very similar GFP [12,33,34]. The average diffusion coefficients are about 3-6 times less in the cellular environment than in water [34].
Previously, nuclear envelope permeabilities have been determined mostly by fitting an analytical solution to a measured average fluorescence intensity in the nucleus after photobleaching [35,36]. The solution assumes that diffusion in the nucleus and cytoplasm is fast compared to diffusion through the nuclear envelope, such that the concentrations  in the nucleus and cytoplasm are approximately uniform. A previously reported value of the permeability for enhanced green fluorescent protein EGFP, determined in this manner, is 0.011 μm/s in COS7 cells [36]. We obtained a clearly higher permeability value and some reasons for the discrepancy are discussed below.
Besides the fact that the permeability of the nuclear envelope is probably somewhat different between two cell lines, an additional explanation might be that the above approximation used in Refs. [35,36] exaggerates the flux through the nuclear envelope (since it assumes a higher concentration difference over the membrane than there actually is). This leads to an underestimation of the nuclear envelope permeability that is obtained by fitting the model to an experiment. Another assumption used in [36] was that the cytoplasm concentration stayed constant, i.e., the cytoplasm volume was large compared to that of the nucleus. This approximation also further underestimates the permeability.
Regions of the cell where these approximations prove especially problematic (for cells on a hard substrate) are the bottom and the top of the nucleus. There the nucleus is usually very close to the plasma membrane of the cell [see Fig. 7

(d)]
and, because of the restricted space, it is more difficult for the cell to replace permeated molecules to retain an efficient concentration difference between the two compartments.

VIII. CONCLUSION
We reported here an analysis of modeling diffusion through semipermeable membranes. We intentionally avoided assigning intrinsic properties such as the diffusion coefficient and particle concentration to the membrane itself, but rather treated the membrane as a semipermeable boundary between two regions. In particular, the transmission boundary conditions were imposed at the membrane. This approach is especially useful in cases where the membrane is thin compared to the overall size of the system. To facilitate a numerical modeling of diffusion through thin membranes, we then presented various discrete treatments of the transmission boundary conditions. First we presented a direct discretization of the conditions at the macroscopic level of description. A mesoscopic discretization, the partial-bounceback scheme, was then shown to be a consistent, alternative approach for enforcing the transmission boundary conditions. Furthermore, an expression for the permeability derived for a mesoscopic LB scheme, was shown to be in accordance with a microscopic, single-particle stochastic diffusion model. Hence, the mesoscopic parameter φ provided a bridge between the microscopic and macroscopic descriptions, and could, for example, be used as a coupling parameter between microscopic and macroscopic solvers in a multiscale simulation effort. Moreover, the proposed lattice-Boltzmann scheme for membrane diffusion is directly applicable to advection-diffusion problems as well as to a hydrodynamic transport across a membrane.
To demonstrate the applicability of the membrane scheme, we simulated nuclear transport in a real cellular geometry imaged by fluorescence microscopy. In this application of the LB membrane scheme, together with a prescribed computational procedure, we numerically determined a value for the permeability of the nuclear envelope to EYFP. Since the number of methods to study nuclear envelope permeability has been limited, our method has the potential to be an important addition to the existing tools.

ACKNOWLEDGMENT
This work was financed by the Jane and Aatos Erkko foundation.

APPENDIX A: ANALYTICAL SOLUTION TO THE DIFFUSION PROBLEM
Consider a diffusion equation, and transmission conditions and initial condition gives us a dimensionless form of the problem Eqs. (A1)-(A4), where γ = LP /D. This problem can be solved by separation of variables, ρ(ξ,τ ) = v(ξ ) w(τ ). For function v we find an eigenvalue problem, and w satisfies the equation,
The transmission condition in Eq. (A5) leads to an eigenvalue equation, Eigenvalues μ 2k are given by while μ 2k+1 are given by the solutions of the equation, The corresponding eigenfunctions are given by From general Sturm-Liouville theory we know that these eigenfunctions are orthogonal with respect to the L 2 inner product, and they form a complete (orthogonal) basis [37]. The solution of the problem of Eq. (A6) can be expressed in the form, × sgn(ξ ) cos(μ 2k+1 ξ ) + 2γ μ 2k+1 sin(μ 2k+1 ξ ) .

APPENDIX B: LATTICE-BOLTZMANN SCHEME OF THE PARTIAL-BOUNCEBACK METHOD
Pure diffusion processes, involving isotropic diffusion, can be simulated, e.g., by using the standard lattice-Boltzmann evolution equation with the BGK-collision operator [11]: where τ is a relaxation time, ρ denotes the local density at a given instant, and suitable weight coefficients for the equilibrium function, f eq i = w i ρ, are given in Table I. The lattice spacing δr, and discrete time step δt, define a reference  The diffusion coefficient can be derived, e.g., by the Chapman-Enskog procedure: Perhaps the simplest model of diffusion through a semipermeable membrane is a one-dimensional random walk in discrete space and time with a special site for the membrane. Such a walker takes steps to the left and right with equal probability except at the special site, where it either goes through the membrane (with probability φ/2,) or is reflected back with probability 1 − φ/2). Hence the membrane vanishes in the limit φ → 1 (i.e., fraction of pores → 1). If N t denotes the number of collisions with the membrane up to time t, the first-passage time through the membrane obeys so the number of visits to the membrane site exceeds a geometric random variable with parameter φ/2 at the time of first passage through the membrane. In other words, the mean number of visits N T to the special site up to and including the first passage equals 2/φ. In particular, two visits on the average are needed to cross a completely permeable membrane (limit φ → 1), and the required number of visits diverges proportional to the first inverse power of φ as φ → 0.
Next we consider a microscopic model continuous in space and time, the probability density of which satisfies the macroscopic diffusion equation Eq. (A1) with the flux boundary condition, Eq. (A3). Following preprint [23], we define the (one-dimensional) snapping-out Brownian motion X t to be a process on the real line that obeys the usual Brownian dynamics except at the origin where the membrane 043309-13 is located. Passing through the membrane consists of a single rate-limiting step given that the particle is performing a search at the immediate vicinity of the membrane (this corresponds to a reaction distance in chemical-reaction applications). The amount of time spent performing the search is measured by the local time, at the membrane, and hence the probability that the firstpassage time T through the membrane is greater than t (i.e., the particle has not passed through the membrane before that) equals where τ c = lim w→0 w/D(w) is the characteristic time scale that depends on the membrane (w is the membrane width). We remark that finding τ c directly from its definition as a limit would require a model for the microscopic structure of the interface. Instead, we now relate time scale τ c to the pore fraction φ by analogy to the random walk. In particular, the number of visits to the membrane site can be considered as the local time of the random walk at the membrane (see Ref. [38] for a convergence proof), and therefore τ c should diverge as ∼1/φ in the limit φ → 0. On the other hand, the first passage can occur immediately upon first contact with the membrane (the first time that the local time L t > 0). In particular, this must hold true in the limit of a completely permeable membrane, so we can assume that τ c → 0 as φ → 1. The only sensible function of φ satisfying these requirements is τ c = 1/φ − 1. Hence the microscopic permeability, defined as the inverse of the mean first passage time, coincides with the macroscopic permeability provided by the Chapman-Enskog analysis: