key: cord-0735958-j976f476
authors: Bressloff, Paul C.
title: Local accumulation time for diffusion in cells with gap junction coupling
date: 2022-01-05
journal: Phys Rev E
DOI: 10.1103/physreve.105.034404
sha: aedd235862aac05b2d351b7e9e9594b2a296e943
doc_id: 735958
cord_uid: j976f476

In this paper we analyze the relaxation to steady-state of intracellular diffusion in a pair of cells with gap-junction coupling. Gap junctions are prevalent in most animal organs and tissues, providing a direct diffusion pathway for both electrical and chemical communication between cells. Most analytical models of gap junctions focus on the steady-state diffusive flux and the associated effective diffusivity. Here we investigate the relaxation to steady state in terms of the so-called local accumulation time. The latter is commonly used to estimate the time to form a protein concentration gradient during morphogenesis. The basic idea is to treat the fractional deviation from the steady-state concentration as a cumulative distribution for the local accumulation time. One of the useful features of the local accumulation time is that it takes into account the fact that different spatial regions can relax at different rates. We consider both static and dynamic gap junction models. The former treats the gap junction as a resistive channel with effective permeability $mu$, whereas the latter represents the gap junction as a stochastic gate that randomly switches between an open and closed state. The local accumulation time is calculated by solving the diffusion equation in Laplace space and then taking the small-$s$ limit. We show that the accumulation time is a monotonically increasing function of spatial position, with a jump discontinuity at the gap junction. This discontinuity vanishes in the limit $mu rightarrow infty$ for a static junction and $beta rightarrow 0$ for a stochastically-gated junction, where $beta$ is the rate at which the gate closes. Finally, our results are generalized to the case of a linear array of cells with nearest neighbor gap junction coupling.

Gap junctions are small nonselective channels that provide a direct diffusion pathway between neighboring cells. They are formed by the head-to-head connection of two hemichannels or connexons, one from each of the two coupled cells [1, 3, 2] . Gap junctions are prevalent in most animal organs and tissues, providing a mechanism for both electrical and chemical communication between cells. Electrical coupling is particularly important in cardiac muscle, where the efficient transmission of electrical signals allows the heart muscle cells to contract in unison. Gap junctions or electrical synapses are also found throughout the central nervous system [4] . Direct chemical communication between cells occurs through the transmission of small second messengers, such as inositol triphosphate (IP3) and calcium (Ca 2+ ). An important example of long-range chemical signaling mediated by gap junctions is the propagation of intercellular Ca 2+ waves [5] .

One of the characteristic properties of a gap junction is its effective channel permeability µ. Mathematically speaking, the high resistance to diffusive flow generates a jump discontinuity ∆u of molecular concentration across the gap junction such that the diffusive flux through the channel is given by J = −µ∆u. A number of studies have analyzed the diffusion equation for a one-dimensional (1D) array of cells with nearest-neighbor coupling [6, 7, 8] . In particular, one can calculate the steady-state concentration in each cell by assuming that there is a constant diffusive flux J 0 through the cellular array; the flux is then determined self-consistently by solving the resulting boundary value problem. The latter includes a set of interior boundary conditions that combine jump discontinuities in the concentration at the gap junctions with flux continuity conditions. Analogous to the opening and closing of ion channels [9] , gap junctions can be both voltage-gated and chemically-gated [10, 11] . This has motivated a stochastic model in which each gap junction is treated as a randomly fluctuating gate that switches between an open and a closed state [12] . Solving the resulting firstorder moment equations for the stochastic concentrations and fluxes in steady state generates an effective channel permeability that depends on the kinetics of the stochastic gate, even though the channel is fully conducting when in the open state.

One advantage of analytical models is that they provide explicit expressions for the steady-state flux J 0 . The latter can be used to extract the effective diffusivity D e for the intracellular transfer of molecules via gap junctions, which can then be compared with experimental data. However, D e is a lumped parameter that depends on both the cytoplasmic diffusivity D and the junctional permeability µ. In order to separate these two biophysical parameters, it is necessary to include additional information about the diffusion process, such as the time-dependent approach to steady-state [7, 8] .

In this paper, we extend previous studies of static and dynamic gap junctions by investigating the relaxation to steady state in terms of the so-called local accumulation time. The latter is commonly used to estimate the time to form a protein concentration gradient during morphogenesis [13, 14, 15, 16, 17] . The basic idea is to treat the fractional deviation from the steady-state concentration as a cumulative distribution for the local accumulation time. One of the useful features of the local accumulation time is that it takes into account the fact that different spatial regions can relax at different rates. (This contrasts with a global measure of the relaxation rate based on the principal nonzero eigenvalue of the negative Laplacian.) In addition, for linear diffusion problems, the mean accumulation time can be calculated by solving the diffusion equation in Laplace space and then taking the small-s limit. This avoids the difficulty in obtaining the full time-dependent solution by evaluating the inverse Laplace transform [7] .

The structure of the paper is as follows. In section II we consider a pair of cells coupled by a single static gap junction and give a definition of the local accumulation time. We solve the resulting diffusion equation in Laplace space, and then use this to calculate both the steady-state flux and the corresponding accumulation time. We show that the latter is a monotonically increasing function of spatial position, with a jump discontinuity at the gap junction. This discontinuity vanishes in the limit µ → ∞. In section III we consider the more complicated example of a pair of cells connected via a stochastically-gated gap junction. Defining the steady-state flux and accumulation time in terms of the solution to the first-order moment equations, we obtain analogous results to the static case. Finally, in section IV our analysis is generalized to the case of a linear array of cells with nearest neighbor gap junction coupling. Throughout the paper we fix the length-scale by taking the size of a cell to be L = 1. The time-scale is then given by L 2 /D.

We begin by considering a pair of identical cells coupled by a single gap junction, as shown in Fig. 1 . Following previous studies [7, 8, 12] , we treat each cell as a one-dimensional compartment of length L and represent the gap junction as a resistive pore with some permeability µ. For the moment we assume that µ is given; a mechanism for generating µ based on stochastic gating will be considered in section 3; see also Ref. [12] . Let u j (x, t), 0 < x < L, denote the concentration of diffusing particles in the j-th cell, j = 1, 2. We then have a pair of diffusion equations

These are supplemented by the interior boundary conditions

where µ is the effective channel permeability, and the exterior boundary conditions

The interior boundary conditions ensure continuity of flux across the gap junction, which depends on the difference in concentrations on either side of the gap junction. A constant concentration difference is maintained between the exterior boundaries of the two cells. (One way to motivate a 1D model is to assume that the external concentrations are uniform with respect to the other spatial dimensions of the cells, and that there is an array of uniformly distributed gap junctions at the interior boundary that can be reduced to a single effective gap junction [8] .) It is straightforward to determine the steady-state concentrations u * j (x), j = 1, 2, by setting all time derivatives to zero and assuming that there is a constant diffusive flux J 0 . The flux is then determined self-consistently by solving the resulting boundary value problem [8] . In particular, one finds that

where D e is an effective diffusion coefficient:

.

(2.5)

As highlighted elsewhere [6, 7] , the effective diffusivity D e can be compared with experimental measurements of tagged particles such as fluorescent probes. However, experimentalists are typically interested in extracting separate values for µ and D, whereas the permeability and diffusivity are lumped together in the expression for D e . As an alternative, one could determine the full timedependent solution by solving Eq. (2.1) in Laplace space and then evaluating the inverse Laplace transform. The resulting solution could then be fit to experimental data [7] . In this paper, rather than obtaining the full time-dependent solution, we focus on a particular characterization of the relaxation to steady state known as the local accumulation time. The latter has previously been developed within the context of diffusion-based morphogenesis [13, 14, 16, 17] , but has more recently been applied to intracellular protein gradient formation [18] and to diffusion processes with stochastic resetting [19] . The accumulation time is easier to calculate than the full time-dependent solution, and provides a more compact representation of the relaxation to steady state. Consider the initial conditions u j (x, 0) =ū j (x) and define

as the fractional deviation of the concentration in the j-th cell from steady state. Assuming that there is no overshooting, 1 − Z j (x, t) can be interpreted as the fraction of the steady-state concentration that has accumulated at x by time t. It follows that −∂ t Z j (x, t)dt is the fraction accumulated in the interval [t, t + dt]. The accumulation time T j (x) at position x ∈ (0, L) is then defined as [13, 14, 16] :

For more complicated diffusion problems it is more convenient to calculate the accumulation time in Laplace space. Let u(x, s)

and, hence

(2.9)

Laplace transforming Eq. (2.1) and the associated boundary conditions gives

for j = 1, 2, with k = s/D, the initial conditions u j (x, 0) =ū j (x), the interior boundary conditions For simplicity, we assume that the cells are initially empty so thatū j ≡ 0. The general solution within each cell then has the form

Substituting these solutions into the various boundary conditions generates four equations for the four unknown coefficients A 1,2 , B 1,2 :

Expressing the solutions (2.13) in terms of B 2 , we have

The solution in Laplace space can now be used to determine the steady-state concentrations and flux according to Eq. (2.8), noting that k = s/D → 0 as s → 0. We thus find from Eqs. (2.16) and (2.17) that

We also recover the flux Eq. (2.4). Although it is simpler to solve the steadystate equations directly, the advantage of working in Laplace space is that one can determine the accumulation time (and other quantities) without having to solve another boundary value problem.

The formula (2.9) implies that we need to evaluate the first derivative of F j (x, s) = s u j (x, s) in the limit s → 0. It can be seen from Eqs. 2.16) and (2.17) that F j (x, s) only depends on s via its k-dependence. Therefore, for fixed x,

This will be non-singular in the limit s → 0 provided that the Taylor expansion of F j about k = 0 is of the form

This is indeed found to be the case, see appendix A, so that

.

(2.21)

Eqs. (2.16) and (2.17a,b) yield the k-dependent functions

and 

and

In contrast to the steady-state flux, the accumulation time could be used to extract values for both D and µ, for example. Consider the composite accumulation time

where Θ(x) is the Heaviside function and T j (x), j = 1, 2, are given by Eqs.

(2.24) and (2.25), respectively. In Fig. 2 we show example plots of T (x) for different diffusivities D and µ = L = 1. As expected, the accumulation time is a strictly monotonically increasing function of x and a decreasing function of D.

In addition, we see that there is a jump discontinuity in T (x) at the gap junction connecting the two cells. The corresponding results for various permabilities µ are shown in Fig. 3 . In the limit µ → ∞ the gap junction no longer restricts the diffusion of molecules and one simply has a domain of length 2L. Taking the limit µ → ∞ in Eqs. (2.24) and (2.25) yields the following expression for the composite accumulation time:

We now consider the more complicated problem of a a pair of cells connected by a stochastically-gated (dynamic) gap junction as formulated in Ref. [12] , see Assume that transitions between the two states n = 0, 1 are described by the two-state Markov process,

The random opening and closing of the gate means that particles diffuse in a random environment according to the piecewise deterministic equations

with static exterior boundary conditions

and n(t)-dependent boundary conditions on the common interior boundary:

That is, when the gate is open there is continuity of the concentration and the flux across the common boundary, whereas when the gate is closed the righthand boundary of cell 1 and the left-hand boundary of cell 2 are reflecting. For simplicity, we assume that the diffusion coefficient is the same in both compartments so that the piecewise nature of the solution is solely due to the switching gate.

Averaging the concentrations with respect to the gate dynamics, we introduce the decompositions

It can be shown that the components satisfy the first-order moment equations [12] ∂V j,

for x ∈ (0, L) and j = 1, 2, with the exterior boundary conditions

and the interior boundary conditions

The boundary condition (3.7a) follows from noting that if n(t) = n and x = 0, then u 1 (0, t) = η with probability one, and thus

Following along similar lines to the static gate, we analyze the diffusion equation in Laplace space. In particular, Eqs. (3.6a,b) become

for x ∈ (0, L) and j = 1, 2; the boundary conditions are the same. From the interior boundary conditions (3.7c), we set

with Γ(s) to be determined later by imposing V 1,0 (L, s) = V 2,0 (0, s). Adding equations (3.6a) and (3.6b) then gives

with the boundary conditions

The corresponding solutions are

We can now substitute for U j (x, s) in equation (3.9b), say, and determine the components V j,1 (x, s) for the two cells.

Introduce the Green's functions G j (x, y; s) according to

with the boundary conditions

and

The Green's functions have the explicit expressions

with A = ξ cosh ξL. The solutions of (3.9b) can then be expressed in terms of the Green's functions as follows:

Finally, the unknown function Γ(s) is determined by requiring that V 1,0 (L, s) = V 2 (0, s) (continuity of the concentration when the gate is open), that is,

with U j given by Eqs (3.12a,b) and V j,1 given by Eqs. (3.17a,b ).

As a check on the above analysis, we begin by deriving the mean steady-state flux J 0 , which was previously obtained by solving the steady-state first moment equations directly [12] . Eqs. (3.11) imply that

Multiplying Eqs. (3.12a,b) by s and taking the limit s → 0 gives

Next, multiplying Eq. (3.17a) by s and taking the limit s → 0 shows that

Multiplying Eq. (3.13) by y m for j = 1, m = 0, 1, and integrating with respect to y, we find that

with ξ 0 = (α + β)/D. We thus obtain the result It is useful to discuss a few limiting cases. First, in the fast switching limit ξ 0 → ∞, we have J 0 → ηD/2L, µ e → ∞ and Eq. (3.20) reduces to the continuous steady-state solution

The mean flux through the gate is the same as the steady-state flux without a gate. On the other hand, for finite switching rates the mean flux J 0 is reduced. In the limit α → 0 (gate always closed), J 0 → 0 so that U * 1 (x) = η for x ∈ [0, L) and U * 2 (x) = 0 for x ∈ (0, L].

In the case of a dynamic gate we define the local accumulation time in terms of the first-order moments of the concentration. That is, assuming that the initial concentrations are zero, we take

where, see Eq. (3.4),

The corresponding accumulation times are then defined according to

with F j (x, s) = s U j (x, s). Note that it does not make sense to define an accumulation time for a single realization of the stochastic gate, since one cannot define a steady-state density, that is, lim t→∞ u j (x, t) does not exist. Eqs. (3.12a) and (3.12b) imply that

and

In appendix B we show that

for small s, with J 0 given by Eq. (3.24), so that

The calculation of J 1 is presented in appendix B, which yields the expression where

In terms of the effective permeability µ e , we have

, (3.34) and

(3.35)

The accumulation times can now be determined by substituting Eqs. (3.20) and (3.32) into (3.28) with J 0 and J 1 given by Eqs. (3.24) and (3.33), respectively. In Fig. 5 we plot the composite accumulation time defined by Eq. (2.26) as a function of x ∈ [0, 2L] for various values of the closing rate β. In the limit β → 0 (ρ 1 → 0) the gate is always open, and we recover the continuous accumulation time of a static gap junction in the limit µ → ∞, see also Fig. 3 . As β increases from zero, however, the accumulation time develops a discontinuity at the junction between the two cells, analogous to the static gap junction with finite permeability. Finally, in Fig. 6 we show analogous plots for different diffusivities. As expected, increasing D reduces the accummulation time. . We will focus on the static case, since the number of discrete states of N independently gated dynamic junctions grows as 2 N , which considerably complicates the analysis [12] . Suppose that we label the cells by an integer ℓ, ℓ = 1, . . . , N , and take the length of each cell to be L. Let u ℓ (x, t) for x ∈ (0, L) denote the particle concentration within the interior of the ℓ-th cell, and consider the diffusion equation

At each of the intercellular boundaries, the concentration is discontinuous due to the permeability of the gap junctions. The generalization of Eq. 

Laplace transforming Eq. (4.1) and the associated boundary conditions gives Substituting these solutions into the various boundary conditions generates the following hierarchy:

Eqs. (4.8b,c) hold for ℓ = 1, . . . , N − 1.

Following along analogous lines to Ref. [7] , we rewrite Eqs. (4.8b,c) as the matrix equation

Let (λ ± , v ± ) denote the eigenpairs of the non-symmetric matrix M and set

(4.11) (For ease of notation, we drop the explicit dependence on k = s/D.) Substituting the eigenfunction expansion into the matrix equation and iterating shows that

In addition, Eqs. (4.8a,d) imply that

(4.14)

We have introduced the dual vectors v j such that

Therefore, setting ℓ = N in Eq. (4.12) yields a pair of equations for the two remaining unknowns B 1 and B N :

Eliminating the coefficient B N we find after some algebra that B 1 has the solution

with

(Recall that the eigenpairs (λ ± , v ± ) depend on k.) Having found B 1 , all of the coefficients {A ℓ , B ℓ , ℓ = 1, . . . , N } are then determined from Eqs. (4.11), (4.12) and (4.13a). In addition, we can now calculate the steady-state concentrations using

and

where f (x) = (e −kx , e kx ) ⊤ . In particular,

and the corresponding steady-state flux is

It remains to calculate the eigenfunctions and eigenvalues of the matrix M(k) and to determine the function g(k). We find that

(4.24b)

The normalization factors N ± ensure the conditions v ± · v ± = 1. However, they cancel out in Eq. (4.18). One important observation is that

It also follows from Eq. (4.14) that Σ ± (k) = Σ ± (0) + O(k) with

Moreover,

The above analysis implies that g(k) = g −1 /k + O(1) and thus 

We thus recover the classical result for the steady-state flux through the coupled system shown in Fig. 7 [7, 8, 12] :

Introducing the effective diffusion coefficient D e according to 

Deriving an explicit expression for the accumulation time is nontrivial for general N . For the sake of illustration, consider the accumulation time T 1 (x), x ∈ [0, L], in the first cell. In this case,

with g(k) given by Eq. (4.18). In Fig. 8 we plot T 1 (x) as a function of x for various cell numbers N . The limit on the right-hand side of Eq. (4.33) is obtained by substituting for g(k) using Eq. (4.18) and numerically evaluating the the first cell.

In order to determine the corresponding accumulation time T ℓ (x) in the ℓ-th cell, we use Eqs. (4.12), (4.13a) and (4.21) . That is, setting

(4.34)

We have used the result B 1 = ηg(k)/s. Note that λ ± , Λ ± and F ± (x) are also functions of k, which means that s u ℓ (x, s) only depends on s via its dependence on k. The generalization of Eq. (4.33) is thus

Define the composite accumulation time

where χ ℓ (x) = 1 for x ∈ [(ℓ − 1)L, ℓL) and is zero otherwise. In Fig. 9 we plot T (x) for the first five cells in a linear array of size N = 10 and various permeabilities µ. As in the case of two cells, see Fig. 3 , T (x) is a monotonically increasing function of x with jump discontinuities at the gap junctions. In the limit µ → ∞ the discontinuities vanish and T (x) becomes a continuous function of x.

In this paper we calculated the local accumulation time for intracellular diffusion in cells with gap-junction coupling. We considered a pair of cells connected by either a static or a dynamic gap junction. In both cases, we showed that the accumulation time is a monotonically increasing function of spatial position with a jump discontinuity at the gap junction. This discontinuity vanished in the limit µ → ∞ for a static junction with permeability µ and in the limit β → 0 for a stochastically-gated junction with rate of closing β. We also extended our analysis of static gap junctions to the case of a linear array of cells with nearest neighbor gap junction coupling. In contrast to the expressions for the steady-state flux, the accumulation times did not simply depend on the effective permeability and cytoplasmic diffusivity via a lumped parameter given by an effective diffusivity. One limitation of our analysis of a stochastically gated gap junction in section III is that we only considered the first-order moment equations for the mean stochastic concentrations E[u j (x, t)]. Mathematically speaking, one can view diffusion in a randomly switching environment such as a stochastically gated gap junction as an example of a piecewise deterministic partial differential equation (PDE). Previously we have shown how one can analyze such a system by discretizing space and constructing the Chapman-Kolmogorov (CK) equation for the resulting finite-dimensional system [20, 12] . The CK equation can then be used to generate a hierarchy of equations for the r-th order moments of the stochastic concentration, which take the form of r-dimensional parabolic PDEs in the continuum limit. Although the diffusing particles are non-interacting, statistical correlations arise at the population level due to the fact that they all move in the same randomly switching environment. That is,

Another simplification of our analysis was to focus on one-dimensional diffusion models. One natural extension of our work would be to consider higher spatial dimensions and more general geometric configurations of cells. In the case of static gap junctions, Keener and Sneyd [8] have analyzed the steadystate flux for a line of two-dimensional cells with gap-junctional openings in the connecting edges. Using symmetry arguments, they showed how the gap junctions along an edge can be lumped into a single effective junction at the center of each edge, whose relative width characterized the degree of clustering of the gap junctions. In particular, they found that clustered gap junctions lead to a much smaller effective diffusion coefficient D e for given gap junctional permeability µ. In future work it would be interesting to investigate how such clustering affects the corresponding local accumulation times. An alternative generalization of one-dimensional diffusion would be to consider diffusion on a tree-like structure with stochastically-gated nodes [21] . A number of biological systems employ branched tree structures in order to distribute nutrients from a single source to many destinations or to gather nutrients from many sources. Examples include plant roots, river basins, neuronal dendrites, and cardiovascular and tracheal systems. The various Green's function moment can be determined from Eq. (3.13). First, multiplying Eq. (3.13) by (L−y) 2 and y 3 , respectively, for j = 1 and integrating with respect to y, we have Combining all of our results, we find that

(1 − ∂ y G 1 (L, 0; 0)) + ∂ y ∂ s G 1 (L, 0; 0), − ρ 1 η∂ s ∂ y G 1 (L, 0; 0) − ρ 1 η L 2 2D ∂ y G 1 (L, 0; 0) = 0, (B.11)

and 

Gap junctions: structure and function

Gap junctions. Cold Spring Harb Perspect

Plasma membrane channels formed by connexins: their regulation and functions

Electrical synapses in the mammalian brain

Intercellular Ca 2+ waves: mechanisms and function

A model for the diffusion of fluorescent probes in the septate giant axon of earthworm: axoplasmic diffusion and junctional membrane permeability

Exact solution of a model of diffusion in an infinite chain or monlolayer of cells coupled by gap junctions

Mathematical Physiology I: Cellular Physiology

Modeling the stochastic gating of ion channels

Gap junction channel gating

Bukauskas A stochastic four-state model of contingent gating of gap junction channels containing two "fast" gates sensitive to transjunctional voltage

Diffusion in cells with stochastically-gated gap junctions

How long does it take to establish a morphogen gradient?

Formation of morphogen gradients: local accumulation time

Physical interpretation of mean local accumulation time of morphogen gradient formation

Local kinetics of morphogen gradients

Critical time scales for advection-diffusion-reaction processes

Protein concentration gradients and switching diffusions

Accumulation time of stochastic processes with resetting

Moment equations for a piecewise deterministic PDE

Diffusion on a tree with stochasticallygated nodes

Carrying out the Taylor expansions of Eqs. 2.22) and (2.23) up to O(k 3 ) we find thatandIt immediately follows that the Taylor expansions have the general form given by Eq. (2.20), namely,Hence, T j (x) = f j (x)/2Du * j (x) and we obtain Eqs. (2.24) and (2.25).

In this appendix we calculate the O(s) contribution to sΓ(s), which is needed in order to determine the local accumulation times in the case of a dynamic gate. Substituting Eqs.