Stochastic Model for the Internal Transfer Kinetics of Cargo in Two-Compartment Carriers

1. Introduction

Mobile nanocarriers are used to deliver drug molecules [1], nucleic acids [2], imaging agents [3], proteins and peptides [4], vaccines [5], and even small metabolites or signaling molecules [6,7]. They come in the form of liposomes, polymeric nanoparticles, micelles, dendrimers, lipid nanoparticles, inorganic nanoparticles, and other engineered nanostructures designed for targeted and controlled delivery [1,8]. Release of cargo from nanocarriers can occur passively through diffusion or gradual degradation of the carrier matrix, or actively in response to specific triggers such as pH, redox gradients, temperature, electric and magnetic fields, or light [9].

The passive release of cargo from a nanocarrier can occur through diffusion out of the carrier into the ambient medium or through collisions with other carriers [10]. For example, a liposome loaded with drug molecules can lose its content over time by diffusion of cargo into the aqueous environment or by random transfer of cargo into a different liposome or another target object upon collision. The acquisition of new cargo through the same two transfer modes is equally possible. Detailed modeling of these mechanisms plays an important role in understanding drug release kinetics and optimizing carrier design as well as in formulating release systems that combine passive and stimulus-responsive release [11,12,13].

Liposomes have been studied extensively as carrier vehicles for drug molecules [14], with more than a dozen formulations being currently approved by the U.S. Food and Drug Administration and the European Medicines Agency for clinical use in cancer therapy, antifungal treatment, pain management, and vaccine delivery [15,16]. The association of a drug with a liposome depends on its hydrophobicity profile. Hydrophilic drugs such as doxorubicin [17] or gemcitabine [18] localize in the aqueous core, with the bilayer acting as a barrier, whereas poorly soluble hydrophobic or amphipathic drugs such as amphotericin B [19] or temoporfin [20] embed in the lipid membrane. Even highly hydrophobic drugs show partial polarity [21], leading to asymmetric leaflet interactions, similar to cholesterol with its hydrophobic backbone buried in the bilayer and hydroxyl group anchored at the headgroup region [22]. Thus, bilayer-associated drug molecules are expected to reside predominantly in two states, being associated either with the external or internal leaflet of the membrane.

The present work aims at contributing to the development of a stochastic model for the transfer kinetics of cargo among mobile nanocarriers. The general scope is illustrated in Figure 1, which shows a schematic representation of carriers with two compartments, external and internal. Although not being displayed in Figure 1, different types of carriers may be present in a general system. Each compartment contains identical sites that are either empty or occupied by a single cargo item.

Upon the collision of two carriers, cargo residing in the external compartment of one carrier can migrate to the external compartment of another carrier. Cargo can also migrate between the two compartments of a given carrier. Hence, the migration from the internal compartment of a carrier to the internal compartment of another carrier must proceed via the two external compartments of the two carriers. Modeling the stochastic time-evolution of a suitable distribution function is a challenging task. In a preceding study [23], we have cast the problem into a set of rate equations and identified an analytic solution in the limit that each carrier contains only one single compartment, which is always filled with a sufficiently large number of cargo but is never even close to be filled completely. We have referred to that limit as the Gaussian regime at low occupation. In the present study, we consider the presence of an external and internal compartment in each carrier and allow for the transfer of cargo between the two. However, we shall not allow for the transfer of cargo between different carriers. This restriction to intra-carrier transport is a significant simplification that renders our stochastic model equivalent to first-order unimolecular reactions [24,25]. Yet, in contrast to having only one “reaction chamber”, our system consists of an ensemble of carriers and thus of one “reaction chamber” for each subpopulation of carriers with fixed total cargo content.

Studying only the internal transfer kinetics of cargo in carriers is a necessary step toward solving the composite problem of allowing for transfer both between and within carriers. Yet, even when ignoring transfer between carriers, the remaining two-compartment stochastic model is of interest because, as we show below, it allows us to express the rate equations as a Fokker-Planck equation in exactly the same limit—the Gaussian regime at low occupation—as in our preceding study [23]. We solve the Fokker-Planck equation and analyze its predictions for three specific examples, including a comparison with numerical solutions of the rate equations. Our work thus paves the way to address the full problem: a multicomponent ensemble of two-compartment carriers with both internal and collision-driven inter-carrier cargo exchange, which we plan to investigate in a future study.

2. Theoretical Model and Discussion

We describe the cargo distribution by the quantity $y_{i,n}\left(t\right)$, denoting the number of carriers that, at a given time t, contain i cargo items in their external compartment and n cargo items in their internal compartment. Each carrier possesses two compartments with identical binding sites, $m_E$ sites in the external and $m_I$ sites in the internal compartment. Every site can either be vacant or host one single cargo molecule. The total number of carriers

$$\begin{array}{c}\hfill N=\sum _{i=0}^{m_E}\sum _{n=0}^{m_I}{y}_{i,n}\left(t\right)\end{array}$$

(1)

is conserved. Cargo can be exchanged between the external and internal compartments of each individual carrier. Hence, the first moments

$$\begin{array}{c}\hfill M_E\left(t\right)=\sum _{i=0}^{m_E}\sum _{n=0}^{m_I}i{y}_{i,n}\left(t\right),M_I\left(t\right)=\sum _{i=0}^{m_E}\sum _{n=0}^{m_I}n{y}_{i,n}\left(t\right)\end{array}$$

(2)

of $y_{i,n}\left(t\right)$ yield time-dependent total numbers of cargo in the external and internal compartment, respectively. The total number $M_E\left(t\right)+M_I\left(t\right)=M$ of cargo in the external and internal compartments of all carriers is conserved. As discussed in the Introduction, we do not consider transfer of cargo among different carriers, implying that ${\sum }_{i=0}^v{y}_{i,v-i}\left(t\right)$ does not depend on time for any fixed v with $0\le v\le m_E+m_I$ and the assumption that $y_{i,n}\left(t\right)=0$ if $i>m_E$ or $n>m_I$.

2.1. Rate Equations

To set up rate equations for the kinetics of intra-carrier transfer of cargo, we display in Figure 2 three carriers that illustrate the processes contributing to the rate of change of $y_{i,n}$.

The population number $y_{i,n}$ increases for a transfer $(i-1,n+1)\to (i,n)$ and for a transfer $(i+1,n-1)\to (i,n)$. Similarly, the population number $y_{i,n}$ decreases for a transfer $(i,n)\to (i-1,n+1)$ and for a transfer $(i,n)\to (i+1,n-1)$. We associate each transfer process with a rate constant, $K_{EI}$ for the transfer of cargo from the external to the internal compartment, and $K_{IE}$ for the transfer of cargo from the internal to the external compartment. The rate constants $K_{EI}$ and $K_{IE}$ reflect passive transport properties across an energy barrier that separates the internal and external compartments. For lipid membranes, permeabilities and flip-flop rates have been studied extensively, both experimentally [26] and via computer simulations [27]. They are found to depend on a multitude of factors, including steric interactions, hydrogen bond formation, membrane composition and degree of asymmetry, the presence of cholesterol, and pH for ionizable drugs. Note that our use of a fixed, single set of rate constants, $K_{EI}$ and $K_{IE}$, neglects cooperativity of transfer processes [28] and variations due to structural and chemical heterogeneities such as locally different compositions in a multi-component lipid vesicle [29,30]. In our present model, transfer events are independent from each other. The rate of change of $y_{i,n}$ takes the form of a chemical Master equation,

$$\begin{array}{c}\hfill {\displaystyle \frac{d{y}_{i,n}\left(t\right)}{dt}}=K_{EI}\left[g_{m_I}^{i+1\to n-1}{y}_{i+1,n-1}\left(t\right)-g_{m_I}^{i\to n}{y}_{i,n}\left(t\right)\right]+K_{IE}\left[g_{m_E}^{n+1\to i-1}{y}_{i-1,n+1}\left(t\right)-g_{m_E}^{n\to i}{y}_{i,n}\left(t\right)\right],\end{array}$$

(3)

with $y_{-1,n}\left(t\right)=y_{m_E+1,n}\left(t\right)=y_{i,-1}\left(t\right)=y_{i,m_I+1}\left(t\right)=0$ for all $0\le i\le m_E$ and $0\le n\le m_I$. Each term in the rate equations also includes a combinatorial factor [31],

$$g_{m_I}^{i\to n}=m_E{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{m_E-1}{i-1}\right)\left(\genfrac{}{}{0pt}{}{m_I-1}{n}\right)}{\left(\genfrac{}{}{0pt}{}{m_E}{i}\right)\left(\genfrac{}{}{0pt}{}{m_I}{n}\right)}}=i\left(1-{\displaystyle \frac{n}{m_I}}\right),g_{m_E}^{n\to i}=m_I{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{m_I-1}{n-1}\right)\left(\genfrac{}{}{0pt}{}{m_E-1}{i}\right)}{\left(\genfrac{}{}{0pt}{}{m_I}{n}\right)\left(\genfrac{}{}{0pt}{}{m_E}{i}\right)}}=n\left(1-{\displaystyle \frac{i}{m_E}}\right),$$

(4)

that accounts for the random statistical occurrence of all possible cargo distributions in a transfer event from the external to the internal compartment (the left of the two equations) and from the internal to the external compartment (the right of the two equations), given the number of cargo in the external and internal compartments is i and n, respectively. Specifically, the factor $g_{m_I}^{i\to n}$ accounts for the number of states $\left({\displaystyle \genfrac{}{}{0pt}{}{m_E-1}{i-1}}\right)\times \left({\displaystyle \genfrac{}{}{0pt}{}{m_I-1}{n}}\right)$ for which a given cargo item can move from an occupied site in the external compartment to an empty site in the internal compartment divided by the total number of available states $\left({\displaystyle \genfrac{}{}{0pt}{}{m_E}{i}}\right)\times \left({\displaystyle \genfrac{}{}{0pt}{}{m_I}{n}}\right)$, multiplied by the number of available sites in the external compartment. For example, $n=m_I$ implies $g_{m_I}^{i\to n}=0$ indicating that transfer into a completely filled internal compartment is not possible. Reasoning for the factor $g_{m_E}^{n\to i}$ is analogous. Because of the presence of the factors $g_{m_I}^{i\to n}$ and $g_{m_E}^{n\to i}$, Eq. 3 can also be referred to as combinatorial Master equation [32]. Based on the rate equations (Eq. 3), we can compute the first moments as defined in Eq. 2, resulting in

$$\begin{array}{c}\hfill {\displaystyle \frac{d{M}_E\left(t\right)}{dt}}=-{\displaystyle \frac{d{M}_I\left(t\right)}{dt}}=K_{IE}{M}_I\left(t\right)-K_{EI}{M}_E\left(t\right)+\left({\displaystyle \frac{K_{EI}}{m_I}}-{\displaystyle \frac{K_{IE}}{m_E}}\right)\sum _{i=0}^{m_E}\sum _{n=0}^{m_I}in{y}_{i,n}\left(t\right).\end{array}$$

(5)

In the limit $m_E\to \infty$ and $m_I\to \infty$ or if the condition $K_{EI}{m}_E=K_{IE}{m}_I$ is fulfilled, the final term on the right-hand side of Eq. 5 vanishes, and the remaining differential equations for $M_E\left(t\right)$ and $M_I\left(t\right)$ yield the well-known exponential behavior of a first-order unimolecular reaction [33],

$${\displaystyle \frac{M_E\left(t\right)}{N}}={\mu }_E+e^{-t/\tau }({\eta }_E-{\mu }_E),{\displaystyle \frac{M_I\left(t\right)}{N}}={\mu }_I+e^{-t/\tau }({\eta }_I-{\mu }_I),$$

(6)

with the characteristic time $\tau =1/(K_{EI}+K_{IE})$ and

$${\mu }_E={\displaystyle \frac{M}{N}}{\displaystyle \frac{K_{IE}}{K_{EI}+K_{IE}}},{\mu }_I={\displaystyle \frac{M}{N}}{\displaystyle \frac{K_{EI}}{K_{EI}+K_{IE}}}.$$

(7)

Note that ${\mu }_E=M_E^{eq}/N$ with $M_E^{eq}=M_E(t\to \infty )$ and ${\mu }_I=M_I^{eq}/N$ with $M_I^{eq}=M_I(t\to \infty )$ denote the number of cargo per carrier in the external and internal compartments, respectively, at equilibrium. Hence, the total amount of cargo per carrier is ${\mu }_E+{\mu }_I=M/N$. Eqs. 6 satisfy the initial conditions ${\eta }_E=M_E(t=0)/N$ and ${\eta }_I=M_I(t=0)/N$ with ${\eta }_E+{\eta }_I=M/N$. Eqs. 6 and 7 are relevant for the present work because below in Sec. Section 2.2 we will adopt the Gaussian limit, where $m_E\to \infty$ and $m_I\to \infty$.

Let us return to the rate equations (Eq. 3). In equilibrium, for $t\to \infty$, the function

$$y_{i,n}^{eq}=y_{i,n}(t\to \infty )=g(i+n)\left({\displaystyle \genfrac{}{}{0pt}{}{m_E}{i}}\right)p^i{(1-p)}^{m_E-i}\left({\displaystyle \genfrac{}{}{0pt}{}{m_I}{n}}\right)q^n{(1-q)}^{m_I-n}$$

(8)

takes on a binomial distribution in both i and n multiplied by a function $g(i+n)$ that depends only on the sum $i+n$ of the two indices and ensures the normalization in Eq. 1 is satisfied. The exact form of $g(i+n)$ depends on the initial distribution $y_{i,n}(t=0)$ as will be discussed below. The probabilities p and q in Eq. 8 relate to the rate constants $K_{EI}$ and $K_{IE}$, the number of sites $m_E$ and $m_I$, and the overall number of cargo per carrier $M/N$ through the two equations

$$m_E{K}_{EI}p(1-q)=m_I{K}_{IE}(1-p)q,p{m}_E+q{m}_I={\displaystyle \frac{M}{N}}.$$

(9)

The equation on the left, which is a consequence of detailed balance [32], results from inserting $y_{i,n}^{eq}$ into Eq. 3 with $d{y}_{i,n}/dt=0$, and the equation on the right reflects the conservation of cargo because $M_E^{eq}=Np{m}_E$ and $M_I^{eq}=Nq{m}_I$. The two relationships in Eq. 9 result from taking the first moments of Eq. 8. The probabilities p and q can be calculated explicitly from Eqs. 9. The general result is somewhat cumbersome, but if both $m_E$ and $m_I$ grow large, it simply becomes

$$p={\displaystyle \frac{K_{IE}}{K_{EI}+K_{IE}}}{\displaystyle \frac{M}{m_{E}N}},q={\displaystyle \frac{K_{EI}}{K_{EI}+K_{IE}}}{\displaystyle \frac{M}{m_{I}N}}.$$

(10)

That is, Eqs. 10 become valid in the Gaussian limit, where $m_E\to \infty$ and $m_I\to \infty$. Recalling Eqs. 7, we observe that $p{m}_E={\mu }_E$ and $q{m}_I={\mu }_I$.

2.2. Continuum Representation in the Gaussian Limit at low Occupation

General analytic solutions for the rate equations in Eq. 3 are not known to us, not even in the limit of large $m_E$ and $m_I$. To make progress, we shall consider two approximations. The first is to adopt the Gaussian limit, where both $m_E$ and $m_I$ become large, but such that p and q adopt values anywhere in the range $0<p<1$ and $0<q<1$. It is then convenient to adopt a continuum description $y_{i,n}\left(t\right)\to y(x,z,t)$, where the continuous variables x and z denote the number of cargo in the external and internal compartments, respectively. That is, the discrete variables i and n in $y_{i,n}\left(t\right)$ are replaced by the continuous variables x and z in $y(x,z,t)$.

The moments in Eqs. 1 and 2 then read $N={\int }_{-\infty }^{\infty }dx{\int }_{-\infty }^{\infty }dzy(x,z,t)$, $M_E\left(t\right)={\int }_{-\infty }^{\infty }dx{\int }_{-\infty }^{\infty }dzxy(x,z,t)$, and $M_I\left(t\right)={\int }_{-\infty }^{\infty }dx{\int }_{-\infty }^{\infty }dzzy(x,z,t)$. The binomial distribution in Eq. 8, which is adopted in equilibrium, turns into a Gaussian,

$$\begin{array}{c}\hfill y^{eq}(x,z)=g(x+z){\displaystyle \frac{e^{-\frac{{(x-m_{E}p)}^2}{2m_{E}p(1-p)}}}{\sqrt{2\pi m_{E}p(1-p)}}}{\displaystyle \frac{e^{-\frac{{(z-m_{I}q)}^2}{2m_{I}q(1-q)}}}{\sqrt{2\pi m_{I}q(1-q)}}}.\end{array}$$

(11)

Instead of the discrete variable $i+n$ in Eq. 8, we employ the corresponding continuous variable $x+z$ as argument of the function $g(x+z)$. It is convenient to re-express the function

$$g(x+z)={\displaystyle \frac{1}{2}}f(x+z){\displaystyle \frac{\sqrt{2\pi [m_{E}p\left(1-p)+m_{I}q(1-q)\right]}}{e^{-\frac{{(x-m_{E}p+z-m_{I}q)}^2}{2\left[m_{E}p(1-p)+m_{I}q(1-q)\right]}}}}$$

(12)

in terms of a function $f(x+z)$. Eq. 12 can be viewed as the definition of the new function $f(x+z)$, given $g(x+z)$ is known. The advantage of using $f(x+z)$ is the simple normalization

$${\displaystyle \frac{1}{2}}\underset{-\infty }{\overset{\infty }{\int }}dvf\left(v\right)=N$$

(13)

that this function must fulfill in order to ensure Eq. 1 is satisfied.

In the Gaussian limit, Eq. 3 can be subjected to a series expansion of $i+1\to x+\Delta x$ and $n+1\to z+\Delta z$ in terms of $\Delta x$ and $\Delta z$ up to first order followed by setting $\Delta x=\Delta z=1$, leading to

$${\displaystyle \frac{dy}{dt}}=\left({\displaystyle \frac{\partial }{\partial x}}-{\displaystyle \frac{\partial }{\partial z}}\right)\left[-K_{IE}\left(1-{\displaystyle \frac{x}{m_E}}\right)zy+K_{EI}\left(1-{\displaystyle \frac{z}{m_I}}\right)xy\right].$$

(14)

This equation is of first order, thus containing drift but no fluctuation terms. Solutions of Eq. 14 would describe how initial delta-peaks would move in time without widening due to fluctuations. Fluctuations would enter Eq. 14 as terms associated with second-order derivatives. Yet, these terms do not emerge correctly from extending the first-order expansion that leads from Eq. 3 to Eq. 14 to a second order expansion. To identify the fluctuation terms, we adopt a second approximation, the low occupation limit, where the number of occupied sites in the external and internal compartment is always much smaller than the number of available sites. The corresponding conditions $M\ll m_{E}N$ and $M\ll m_{I}N$ imply $p\ll 1$ and $q\ll 1$. Note that ${\mu }_E=p{m}_E$ and ${\mu }_I=q{m}_I$ both remain finite. Mathematically, the small occupation limit uses the Gaussian distribution $e^{-{(x-\mu )}^2/\left(2\mu \right)}/\sqrt{2\pi \mu }$ to approximate a Poisson distribution $e^{-\mu }{\mu }^x/x!$ for large mean value $\mu$ [34]. The equilibrium distribution specified in Eqs. 11 and 12 then reads

$$y^{eq}(x,z)={\displaystyle \frac{1}{2}}f(x+z){\displaystyle \frac{\sqrt{2\pi ({\mu }_E+{\mu }_I)}}{e^{-\frac{{(x-{\mu }_E+z-{\mu }_I)}^2}{2({\mu }_E+{\mu }_I)}}}}{\displaystyle \frac{e^{-\frac{{(x-{\mu }_E)}^2}{2{\mu }_E}}}{\sqrt{2\pi {\mu }_E}}}{\displaystyle \frac{e^{-\frac{{(z-{\mu }_I)}^2}{2{\mu }_I}}}{\sqrt{2\pi {\mu }_I}}},$$

(15)

with the function $f(x+z)$ satisfying the normalization in Eq. 13. In the low occupation limit, the terms $x/m_E$ and $z/m_I$ in Eq. 14 are negligibly small, and it becomes straightforward to identify the missing second-order fluctuation term by comparing the stationary solutions of Eq. 14 with Eq. 15. The resulting equation

$$\tau {\displaystyle \frac{dy}{dt}}=\left({\displaystyle \frac{\partial }{\partial x}}-{\displaystyle \frac{\partial }{\partial z}}\right)\left[{\displaystyle \frac{{\mu }_{I}x-{\mu }_{E}z}{{\mu }_E+{\mu }_I}}y+{\displaystyle \frac{{\mu }_E{\mu }_I}{{\mu }_E+{\mu }_I}}\left({\displaystyle \frac{\partial y}{\partial x}}-{\displaystyle \frac{\partial y}{\partial z}}\right)\right]$$

(16)

is a Fokker-Planck equation [32] for the distribution $y(x,z,t)$, with $\tau =1/(K_{EI}+K_{IE})$. We reiterate that the drift term, which is associated with the first derivatives, follows from Eq. 14 in the limit of large $m_E$ and $m_I$ and employs the definitions of ${\mu }_E$ and ${\mu }_I$ in Eq. 7. The fluctuation term, which accounts for the second order derivatives in Eq. 16, ensures the equilibrium distribution $y^{eq}(x,z)=y(x,z,t\to \infty )$ of Eq. 16 is given by Eq. 15.

2.3. Solution of the Fokker-Planck Equation

In order to solve Eq. 16 we introduce the new independent variables $u=x-z$ and $v=x+z$. This is motivated by the conservation of the total number of cargo, v, in each carrier population $y\left(\right(v+u)/2,(v-u)/2,t)$ as u varies. The function $\tilde{y}(u,v,t)=y((v+u)/2,(v-u)/2,t)$ then satisfies the equation

$$\tau {\displaystyle \frac{d\tilde{y}}{dt}}={\displaystyle \frac{\partial }{\partial u}}\left\{\left[u-v{\displaystyle \frac{{\mu }_E-{\mu }_I}{{\mu }_E+{\mu }_I}}\right]\tilde{y}+{\displaystyle \frac{4{\mu }_E{\mu }_I}{{\mu }_E+{\mu }_I}}{\displaystyle \frac{\partial \tilde{y}}{\partial u}}\right\}.$$

(17)

Because we adopted the Gaussian limit, the new variables u and v vary from $-\infty$ to ∞. Also, note that $dxdz=dudv/2$. Hence, solutions of Eq. 17 must conserve the number of carriers $N=1/2\times {\int }_{-\infty }^{\infty }dv{\int }_{-\infty }^{\infty }du\tilde{y}(u,v,t)$ and the total number of cargo $M=1/2\times {\int }_{-\infty }^{\infty }dvv{\int }_{-\infty }^{\infty }du\tilde{y}(u,v,t)$. Let us specify an initial distribution $\tilde{y}(u,v,t=0)={\tilde{y}}_0(u,v)$. Solutions of Eq. 17,

$$\tilde{y}(u,v,t)=\underset{-\infty }{\overset{\infty }{\int }}d{u}^{\prime }G(u,v,t|u^{\prime }){\tilde{y}}_0(u^{\prime },v),$$

(18)

can be constructed using the Green’s function,

$$G(u,v,t|u^{\prime })={\displaystyle \frac{1}{\sqrt{2\pi \sigma \left(t\right)}}}{e}^{-{\displaystyle \frac{{\left[u-\mu (v,t|u^{\prime })\right]}^2}{2\sigma \left(t\right)}}},$$

(19)

with

$$\begin{array}{c}\hfill \mu (v,t|u^{\prime })=u^{\prime }{e}^{-t/\tau }+v{\displaystyle \frac{{\mu }_E-{\mu }_I}{{\mu }_E+{\mu }_I}}\left(1-e^{-t/\tau }\right),\sigma \left(t\right)={\displaystyle \frac{4{\mu }_E{\mu }_I}{{\mu }_E+{\mu }_I}}\left(1-e^{-2t/\tau }\right).\end{array}$$

(20)

The Green’s function $G(u,v,t|u^{\prime })$, defined through Eqs. 19 and 20, satisfies the Fokker-Planck equation (Eq. 17) and is normalized such that ${\int }_{-\infty }^{\infty }duG(u,v,t|u^{\prime })=1$. It initially produces a delta-peak at position $u^{\prime }$ — that is, the Green’s function satisfies the initial condition $G(u,v,t\to 0|u^{\prime })=\delta (u-u^{\prime })$, where $\delta \left(u\right)$ denotes the Dirac delta function. Hence, Eq. 18 reproduces the initial distribution ${\tilde{y}}_0(u,v)$ in the limit $t\to 0$. In the opposite limit, at $t\to \infty$, the Green’s function

$$G(u,v,t\to \infty |u^{\prime })={\displaystyle \frac{1}{\sqrt{2\pi \frac{4{\mu }_E{\mu }_I}{{\mu }_I+{\mu }_E}}}}{e}^{-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\left(u-v\frac{{\mu }_E-{\mu }_I}{{\mu }_I+{\mu }_E}\right)}^2}{\frac{4{\mu }_E{\mu }_I}{{\mu }_I+{\mu }_E}}}}$$

(21)

is independent of $u^{\prime }$, and Eq. 18 thus yields

$$\tilde{y}(u,v,t\to \infty )=G(u,v,t\to \infty |u^{\prime })\underset{-\infty }{\overset{\infty }{\int }}d{u}^{\prime }{\tilde{y}}_0(u^{\prime },v).$$

(22)

The normalization $N=1/2\times {\int }_{-\infty }^{\infty }dv{\int }_{-\infty }^{\infty }du\tilde{y}(u,v,t)$ together with Eq. 13 implies $f\left(v\right)={\int }_{-\infty }^{\infty }du\tilde{y}(u,v,t)$ is conserved. The integral in Eq. 22 is thus ${\int }_{-\infty }^{\infty }d{u}^{\prime }{\tilde{y}}_0(u^{\prime },v)=f\left(v\right)$. Using this, together with Eq. 21 and the definitions $u=x-z$ and $v=x+z$, renders Eq. 22 identical to Eq. 15, thus reproducing the correct equilibrium distribution.

It is also interesting to calculate moments of our solution $\tilde{y}(u,v,t)$ in Eq. 18. The first moment with respect to v yields

$$\begin{array}{ccc}\hfill \langle v\rangle & =& {\displaystyle \frac{1}{2N}}\underset{-\infty }{\overset{\infty }{\int }}du\underset{-\infty }{\overset{\infty }{\int }}dvv\tilde{y}(u,v,t)={\displaystyle \frac{1}{2N}}\underset{-\infty }{\overset{\infty }{\int }}du\underset{-\infty }{\overset{\infty }{\int }}dvv\underset{-\infty }{\overset{\infty }{\int }}d{u}^{\prime }G(u,v,t|u^{\prime }){\tilde{y}}_0(u^{\prime },v)\hfill \\ & =& {\displaystyle \frac{1}{2N}}\underset{-\infty }{\overset{\infty }{\int }}d{u}^{\prime }\underset{-\infty }{\overset{\infty }{\int }}dvv{\tilde{y}}_0(u^{\prime },v)={\displaystyle \frac{M}{N}},\hfill \end{array}$$

(23)

demonstrating that our solution $\tilde{y}(u,v,t)$ indeed conserves M. It can analogously be concluded that higher moments in v are all conserved, implying that distribution functions $y(x,z,t)$ can change with time in the $x-z$ direction but not in the $x+z$ direction. The first moment with respect to u results in

$$\begin{array}{ccc}\hfill \langle u\rangle & =& {\displaystyle \frac{1}{2N}}\underset{-\infty }{\overset{\infty }{\int }}du\underset{-\infty }{\overset{\infty }{\int }}dvu\tilde{y}(u,v,t)={\displaystyle \frac{1}{2N}}\underset{-\infty }{\overset{\infty }{\int }}du\underset{-\infty }{\overset{\infty }{\int }}dvu\underset{-\infty }{\overset{\infty }{\int }}d{u}^{\prime }G(u,v,t|u^{\prime }){\tilde{y}}_0(u^{\prime },v)\hfill \\ & =& {\mu }_E-{\mu }_I+e^{-t/\tau }\left({\eta }_E-{\eta }_I-{\mu }_E+{\mu }_I\right)={\displaystyle \frac{M_E\left(t\right)}{N}}-{\displaystyle \frac{M_I\left(t\right)}{N}},\hfill \end{array}$$

(24)

thus reproducing $M_E\left(t\right)$ and $M_I\left(t\right)$ according to Eq. 6. To transition from the first to the second line of Eq 24, we have made use of ${\mu }_E+{\mu }_I=M/N$ (see Eq. 7), $\tau =1/(K_{IE}+K_{EI})$, ${\eta }_E=M_E\left(0\right)/N$, ${\eta }_I=M_I\left(0\right)/N$, and ${\int }_{-\infty }^{\infty }duuG(u,v,t|u^{\prime })=\mu (v,t|u^{\prime })$ (see Eq. 20). Calculation of the second moment with respect to u gives rise to

$$\begin{array}{ccc}\hfill \langle u^2\rangle ={\displaystyle \frac{1}{2N}}\underset{-\infty }{\overset{\infty }{\int }}du\underset{-\infty }{\overset{\infty }{\int }}dv{u}^2\tilde{y}(u,v,t)& =& {\displaystyle \frac{4{\mu }_E{\mu }_I}{{\mu }_E+{\mu }_I}}\left(1-e^{-2t/\tau }\right)+e^{-2t/\tau }{\langle u^2\rangle }_0\hfill \\ & +& 2{\displaystyle \frac{{\mu }_E-{\mu }_I}{{\mu }_E+{\mu }_I}}\left(e^{-t/\tau }-e^{-2t/\tau }\right){\langle uv\rangle }_0\hfill \\ & +& {\left({\displaystyle \frac{{\mu }_E-{\mu }_I}{{\mu }_E+{\mu }_I}}\right)}^2{\left(1-e^{-t/\tau }\right)}^2{\langle v^2\rangle }_0,\hfill \end{array}$$

(25)

where we have used ${\int }_{-\infty }^{\infty }du{u}^{2}G(u,v,t|u^{\prime })=\mu {(v,t|u^{\prime })}^2+\sigma \left(t\right)$, and where we have defined the initial value (at $t=0$) of the second moment,

$$\begin{array}{c}\hfill {\langle u^2\rangle }_0={\displaystyle \frac{1}{2N}}\underset{-\infty }{\overset{\infty }{\int }}du\underset{-\infty }{\overset{\infty }{\int }}dv{u}^2{\tilde{y}}_0(u,v),\end{array}$$

(26)

with analogous definitions for ${\langle uv\rangle }_0$ and ${\langle v^2\rangle }_0$. Note $\langle v^2\rangle ={\langle v^2\rangle }_0$ is conserved, as discussed above. Eq. 25 demonstrates how the second moment with respect to u propagates in time from the initial ${\langle u^2\rangle }_0$ to the final value,

$$\begin{array}{c}\hfill {\langle u^2\rangle }_{eq}={\displaystyle \frac{1}{2N}}\underset{-\infty }{\overset{\infty }{\int }}du\underset{-\infty }{\overset{\infty }{\int }}dv{u}^2{\tilde{y}}^{eq}(u,v)={\displaystyle \frac{4{\mu }_E{\mu }_I}{{\mu }_E+{\mu }_I}}+{\left({\displaystyle \frac{{\mu }_E-{\mu }_I}{{\mu }_E+{\mu }_I}}\right)}^2{\langle v^2\rangle }_0,\end{array}$$

(27)

which depends on ${\langle v^2\rangle }_0$. That is, if two initial distributions $y^{eq}(x,z)$ differ in their second moment ${\langle v^2\rangle }_0$, then their equilibrium distributions will also be different, given that ${\mu }_E\ne {\mu }_I$. Based on $\langle u^2\rangle$ in Eq. 25 and $\langle u\rangle$ in Eq. 24, the calculation of the variance $〈{\left(u-\langle u\rangle \right)}^2〉=\langle u^2\rangle -{\langle u\rangle }^2$ is straightforward. As pointed out, the corresponding variance in the v-direction, $〈{\left(v-\langle v\rangle \right)}^2〉=\langle v^2\rangle -{\langle v\rangle }^2={\langle v^2\rangle }_0-{(M/N)}^2$ is independent of time.

2.4. Three Specific Examples

The following three examples illustrate simple cases for the calculation of $y(x,z,t)$ from a given initial distribution $y_0(x,z)$, rate constants $K_{EI}$ and $K_{IE}$, and overall cargo-to-carrier ratio $M/N={\mu }_E+{\mu }_I$.

2.4.1. First Example

Assume all N carriers contain exactly the same number of cargo, $M/N$, in each carrier, and all of that is initially contained exclusively in the external compartment, $M_E\left(0\right)=M$ and $M_I\left(0\right)=0$. Our assumptions amount to the initial distribution function ${\tilde{y}}_0(u,v)=2N\delta (u-M/N)\delta (v-M/N)$. The general solution in Eq. 18 then reads $\tilde{y}(u,v,t)=2NG(u,v,t|v)\delta (v-M/N)$. Because of the presence of the delta-peak $\delta (v-M/N)$, it is convenient and sufficient to focus on

$$\begin{array}{c}\hfill \underset{-\infty }{\overset{\infty }{\int }}dv\tilde{y}(u,v,t)=2NG\left(u,{\displaystyle \frac{M}{N}},t|{\displaystyle \frac{M}{N}}\right)={\displaystyle \frac{2N}{\sqrt{2\pi \sigma \left(t\right)}}}{e}^{-{\displaystyle \frac{{\left[u-{\mu }_E+{\mu }_I\left(1-2e^{-t/\tau }\right)\right]}^2}{2\sigma \left(t\right)}}},\end{array}$$

(28)

with $\sigma \left(t\right)$ specified in Eq. 20. Eq. 28 describes how the distribution shifts from the initial delta-peak at $u=M/N$ and $v=M/N$ to a Gaussian at $u={\mu }_E-{\mu }_I$ and $v=M/N$. Expressed in terms of x and z, this shift occurs from the initial $x=M/N$ and $z=0$ to the final $x={\mu }_E=M_E^{eq}/N$ and $z={\mu }_I=M_I^{eq}/N$. Calculation of relevant moments gives rise to $\langle v\rangle =M/N$, $〈{\left(v-\langle v\rangle \right)}^2〉=0$, $\langle u\rangle ={\mu }_E-{\mu }_I\left(1-2e^{-t/\tau }\right)$, and $〈{\left(u-\langle u\rangle \right)}^2〉=\sigma \left(t\right)$. The latter reflects the time evolution of the Green’s function, consistent with our initial distribution ${\tilde{y}}_0(u,v)$ consisting of a delta-peak.

2.4.2. Second Example

We assume an initially Gaussian distribution

$$\begin{array}{c}\hfill y_0(x,z)={\displaystyle \frac{N}{2\pi \sqrt{{\eta }_E{\eta }_I}}}{e}^{-{\displaystyle \frac{{(x-{\eta }_E)}^2}{2{\eta }_E}}}{e}^{-{\displaystyle \frac{{(z-{\eta }_I)}^2}{2{\eta }_I}}},\end{array}$$

(29)

where we recall from Eq. 6 the initial numbers of cargo per carrier, ${\eta }_E=M_E(t=0)/N$ and ${\eta }_I=M_I(t=0)/N$ in the two compartments, with ${\eta }_E+{\eta }_I=M/N$ being the combined number. Using ${\tilde{y}}_0(u,v)=y_0((v+u)/2,(v-u)/2)$, we can calculate the function

$$\begin{array}{c}\hfill f\left(v\right)=\underset{-\infty }{\overset{\infty }{\int }}du{\tilde{y}}_0(u,v)={\displaystyle \frac{2N}{\sqrt{2\pi M/N}}}{e}^{-{\displaystyle \frac{{(v-M/N)}^2}{2M/N}}}.\end{array}$$

(30)

Inserting $f\left(v\right)=f(x+y)$ into Eq. 15 yields the equilibrium distribution

$$\begin{array}{c}\hfill y^{eq}(x,z)={\displaystyle \frac{N}{2\pi \sqrt{{\mu }_E{\mu }_I}}}{e}^{-{\displaystyle \frac{{(x-{\mu }_E)}^2}{2{\mu }_E}}}{e}^{-{\displaystyle \frac{{(z-{\mu }_I)}^2}{2{\mu }_I}}}.\end{array}$$

(31)

Calculating the full kinetic behavior $y(x,z,t)$ from $\tilde{y}(u,v,t)$ according to Eqs. 17-20 leads to somewhat cumbersome expressions but can, of course, be carried out numerically. Here, we focus on the specific case ${\mu }_E={\eta }_I=2/3\times M/N$ and ${\mu }_I={\eta }_E=1/3\times M/N$, where the initial distribution of the cargo in the two compartments is opposite of the preferred one. This leads to the simple analytic result,

$$\begin{array}{ccc}\hfill y(x,z,t)& =& y^{eq}(x,z)e^{-{\displaystyle \frac{N}{4M}}(x+z)\left[e^{-2t/\tau }(x+z)+2e^{-t/\tau }(x-2z)\right]},\hfill \end{array}$$

(32)

where the limiting values for $t=0$ and $t\to \infty$,

$$\begin{array}{ccc}\hfill y_0(x,z)& =& {\displaystyle \frac{N}{2\pi \sqrt{\frac{M}{3N}\frac{2M}{3N}}}}{e}^{-{\displaystyle \frac{3}{2}}{\displaystyle \frac{N}{M}}{\left(x-{\displaystyle \frac{1}{3}}{\displaystyle \frac{M}{N}}\right)}^2}{e}^{-{\displaystyle \frac{3}{4}}{\displaystyle \frac{N}{M}}{\left(z-{\displaystyle \frac{2}{3}}{\displaystyle \frac{M}{N}}\right)}^2},\hfill \\ \hfill y^{eq}(x,z)& =& {\displaystyle \frac{N}{2\pi \sqrt{\frac{M}{3N}\frac{2M}{3N}}}}{e}^{-{\displaystyle \frac{3}{4}}{\displaystyle \frac{N}{M}}{\left(x-{\displaystyle \frac{2}{3}}{\displaystyle \frac{M}{N}}\right)}^2}{e}^{-{\displaystyle \frac{3}{2}}{\displaystyle \frac{N}{M}}{\left(z-{\displaystyle \frac{1}{3}}{\displaystyle \frac{M}{N}}\right)}^2},\hfill \end{array}$$

(33)

agree with Eq. 29 and 31. The left diagram of Figure 3 shows a contour plot of $y(x,z,t)$ according to Eq. 32 for $M/N=100$ at the five different times specified in the legend.

The dashed line displays the path

$$\begin{array}{c}\hfill {\displaystyle \frac{M_E\left(t\right)}{N}}={\displaystyle \frac{M}{3N}}\left(2-e^{-t/\tau }\right),{\displaystyle \frac{M_I\left(t\right)}{N}}={\displaystyle \frac{M}{3N}}\left(1+e^{-t/\tau }\right),\end{array}$$

(34)

along which the maximum of the Gaussian distribution $y(x,z,t)$ moves over time until reaching equilibrium. Eqs. 34 are identical to Eqs. 6, with the specific choices ${\mu }_E=2{\mu }_I=2{\eta }_E={\eta }_I=2M/\left(3N\right)$ of our current example. Because of $M_E\left(t\right)+M_I\left(t\right)=M$, this path is linear, with a slope $-1$ of the dashed line in Figure 3. The probability distributions along all lines of slope $-1$ correspond to a total number of cargo $x+z$. That is, only $u=x-z$ but not $v=x+z$ changes along a line of slope $-1$. Note also $M_E\left(t\right)=M_I\left(t\right)$ in Eq. 34 occurs after a time $t=\tau ln2$. Hence, for the green curves in Figure 3 the cargo is distributed equally among external and internal compartment, $M_E/M=M_I/M=1/2$. The sequence of the five selected times $t/\tau =0$, $ln(4/3)$, $ln\left(2\right)$, $ln\left(4\right)$, and ∞, in Figure 3 corresponds to $M_E/M=4/12$, $5/12$, $6/12$, $7/12$, and $8/12$.

Moments for our specific example can be calculated using the general formalism in Eqs. 23-25 or, equivalently, directly from the solution in Eq. 32, yielding for the variances

$$〈{\left(v-\langle v\rangle \right)}^2〉={\displaystyle \frac{M}{N}},〈{\left(u-\langle u\rangle \right)}^2〉={\displaystyle \frac{M}{N}}\left[1-{\displaystyle \frac{4}{9}}\left(e^{-t/\tau }-e^{-2t/\tau }\right)\right].$$

(35)

As pointed out above, averaging over v yields time-independent results. The corresponding variance in the u-direction is $M/N$ for $t=0$ and $t\to \infty$, but adopts a minimum $8/9\times M/N$ at $t/\tau =ln2$. Hence, the variances for the distributions at $t=0$ (purple) and $t\to \infty$ (red) in the left diagram of Figure 3 are the same in the $x-z$-direction and in the $x+z$-direction. At all intermediate times, the variance in the $x-z$-direction is smaller than that in the $x+z$-direction, with the largest difference for $t=\tau ln2$ (shown in green).

The right diagram in Figure 3 shows cross-sections of $y(x,z,t)$ along the direction of the dashed line in the left diagram, where $x+z=M/N=100$ is fixed. Curves are displayed for the same set of times as in the left diagram, with the same color code. The solid lines correspond to $y(x,z,t)$ in Eq. 32. The open circles show numerical solutions of the corresponding discrete rate equations, Eq. 3, with $K_{EI}=1/\left(3\tau \right)$ and $K_{IE}=2/\left(3\tau \right)$, $N=100$, $M=10,000$, and an initial distribution $y_{i,n}(t=0)=N\left({\displaystyle \genfrac{}{}{0pt}{}{m_E}{i}}\right)p_0^i{(1-p_0)}^{m_E-i}\left({\displaystyle \genfrac{}{}{0pt}{}{m_I}{n}}\right)q_0^n{(1-q_0)}^{m_I-n}$, where $p_0=M/\left(3N{m}_E\right)$ and $q_0=2M/\left(3N{m}_I\right)$. Optimal agreement is expected in the limit $m_E\to \infty$ and $m_I\to \infty$. Our calculations were performed for $m_E=m_I=2500$, which also leads to excellent matching. The comparison in the right diagram of Figure 3 of corresponding numerical solutions for the discrete rate equations in Eq. 3 and analytic solutions of the Fokker-Planck equation in Eq. 16 reinforces the equivalence of the two approaches.

2.4.3. Third Example

Our final example illustrates the addition of another Gaussian distribution to the initial distribution $y_0(x,z)$. That is, we replace the single Gaussian distribution in Eq. 29 by the sum of two Gaussians,

$$\begin{array}{c}\hfill y_0(x,z)={\displaystyle \frac{\alpha N}{2\pi \sqrt{{\eta }_E^{\left(1\right)}{\eta }_I^{\left(1\right)}}}}{e}^{-{\displaystyle \frac{{\left(x-{\eta }_E^{\left(1\right)}\right)}^2}{2{\eta }_E^{\left(1\right)}}}}{e}^{-{\displaystyle \frac{{\left(z-{\eta }_I^{\left(1\right)}\right)}^2}{2{\eta }_I^{\left(1\right)}}}}+{\displaystyle \frac{(1-\alpha )N}{2\pi \sqrt{{\eta }_E^{\left(2\right)}{\eta }_I^{\left(2\right)}}}}{e}^{-{\displaystyle \frac{{\left(x-{\eta }_E^{\left(2\right)}\right)}^2}{2{\eta }_E^{\left(2\right)}}}}{e}^{-{\displaystyle \frac{{\left(z-{\eta }_I^{\left(2\right)}\right)}^2}{2{\eta }_I^{\left(2\right)}}}}.\end{array}$$

(36)

The two diagrams of Figure 4 show $y(x,z,t)$ for two specific choices for the parameters ${\eta }_E^{\left(1\right)}$, ${\eta }_E^{\left(2\right)}$, ${\eta }_I^{\left(1\right)}$, ${\eta }_I^{\left(2\right)}$, and $\alpha$. The rate constants $K_{EI}$ and $K_{IE}$ as well as $M/N$ are the same as in Figure 3.

In the left diagram, the two Gaussians of Figure 3 for $t=0$ and $t/\tau =0.4$ are added, with $\alpha =1/2$ giving each contribution the same weight. That is, ${\eta }_E^{\left(1\right)}=1/3\times M/N$, ${\eta }_I^{\left(1\right)}=2/3\times M/N$, ${\eta }_E^{\left(2\right)}/{\eta }_E^{\left(1\right)}=1+e^{-4/10}$, and ${\eta }_I^{\left(2\right)}/{\eta }_I^{\left(1\right)}=1-e^{-4/10}$ were specified. Contour lines of the resulting initial distribution $y_0(x,z)$ are displayed in purple. The cumulative probabilities $f(x+z)=f\left(v\right)={\int }_{-\infty }^{\infty }du{y}_0(u,v)$ as well as ${\mu }_E$ and ${\mu }_I$ are the same in Figure 3 and the left diagram of Figure 4. Hence, the equilibrium distributions $y^{eq}(x,z)$ (with contour lines shown in red) are the same.

In the right diagram, we separate the two Gaussians along the x-axis by choosing ${\eta }_E^{\left(1\right)}=(1/3-1/10)\times M/N$, ${\eta }_E^{\left(2\right)}=(1/3+1/10)\times M/N$, ${\eta }_I^{\left(1\right)}={\eta }_I^{\left(2\right)}=2/3\times M/N$ as well as $\alpha =1/2$. The purple contour lines represent $y_0(x,z)$. Here, the initial distribution evolves toward a different final distribution (with contour lines shown in red) because $f(x+z)$ for the right diagram of Figure 4 is different from $f(x+z)$ for Figure 3 and the left diagram of Figure 4.

3. Conclusions

This work presents a stochastic model for the kinetics of cargo transfer among two compartments that each contain identical sites but different affinities. We show that rate equations give rise to a Fokker-Planck equation in the continuum and low occupation limits, where the Gaussian and Poisson regimes of a binomial distribution overlap. That implies a large number of cargo must be present yet without ever approaching full occupation of either compartment. Analytic solutions of the Fokker-Planck equation are presented, based on identifying a Green’s function, and discussed in terms of moments and three specific examples. We emphasize that our model represents a step towards rationalizing the kinetics of collision-dominated cargo transfer among mobile two-compartment nanocarriers such as drug molecules associated with the inner or outer leaflets of lipid vesicles. While a preceding study has focused on the inter-carrier transfer of cargo [23], the present one addresses the intra-carrier transport. Future work will combine the two models into one comprehensive approach.