We summarize here the main passages of the generalization of the Fick's first law. We refer the reader to Lecca et al.[25] for a more comprehensive description of the mathematical structures and passages.
A generalization of Fick's law for modeling diffusion
In a chemical system the driving force for diffusion of each species is the gradient of chemical potential µ of this species. The chemical potential of any particular chemical species i is
(1)
where is the standard chemical potential of the species i (i.e. the Gibbs energy of 1 mol of species i at a pressure of 1 bar), R = 8.314 J · K-1 · mol-1 is the ideal gas constant, and T the absolute temperature.
The quantity a
i
is called chemical activity of component i, and it is given by
where γ
i
is the activity coefficient, and c0 a reference concentration, which, for example, could be set equal to the initial concentration. The activity coefficients express the deviation of a solution from ideal thermodynamic behavior and, in general, they may depend on the concentration of all the solutes in the system. For an ideal solution, the limit of γ
i
, which is recovered experimentally at high dilutions is γ
i
= 1. If the concentration of species i varies from point to point in space, then so does the chemical potential. For simplicity, here the case in which there is only a chemical potential gradient in the x direction is taken into account. Chemical potential is the free energy per mole of substance, free energy is the negative of the work, W, which a system can perform, and work is connected to force F acting on the molecules by dW = Fdx. Therefore an inhomogeneous chemical potential is related to a virtual force per molecule of
(3)
where N
A
= 6.022 × 1023 mol-1 is Avogadro's number, k
B
= 1.381 × 10-23 J · K-1 is the Boltzmann constant, and the sum is taken over all species in the system other than the solvent. This force is balanced by the drag force experienced by the solute (F
drag
,
i
) as it moves through the solvent. Drag forces are proportional to speed. If the speed of the solute is not too high in such a way that the solvent does not exhibit turbulence, the drag force can be written as follows
(4)
where f
i
α c
i
is the frictional coefficient, and v
i
is the mean drift speed.
Moreover, if the solvent is not turbulent, the flux, defined as the number of moles of solute which pass through a small surface per unit time per unit area, can be approximated as in the following
i.e. the number of molecules per unit volume multiplied by the linear distance traveled per unit time.
Since the virtual force on the solute is balanced by the drag force (i. e. F
drag,i
= -F
i
), the following expression for the mean drift velocity is obtained
so that Eq. (5) becomes
(6)
where
(7)
are the diffusion coefficients. The Eq. (7) states that, in general, the flux of one species depends on the gradients of all the others, and not only on its own gradient. However, here it is supposed that the chemical activity a
i
depends only weakly on the concentrations of the other solutes, i.e. it is assumed that D
ij
≈ 0 for i ≠ j and that Fick's laws still holds [26]. Let D
i
denote D
ii
. It is still generally the case that D
i
depends on c
i
in sufficiently concentrated solutions since γ
i
(and thus a
i
) has a non trivial dependence on c
i
[26]. There is only one very special case, namely that of an ideal solution with γ
i
= 1, in which the diffusion coefficient, D
i
= k
B
T/f
i
, is constant. In order to find an analytic expression for the diffusion coefficient, D
i
, in terms of the concentration, c
i
, let us consider that the rate of change of concentration of the substance i due to diffusion is given by
Substituting Eq. (7) into Eq. (6), and then substituting the obtained expression for J
i
into Eq. (8), give
so that
Let c
i
,
k
denote the concentration of a substance i at coordinate x
k
, and l = x
k
- x
k
-1 the distance between adjacent mesh points. The derivative of c
i
with respect to x calculated at is
(9)
By using Eq. (9) into Eq. (6) the diffusive flux of species i midway between the mesh points, is obtained:
(10)
where is the diffusion coefficient midway between the mesh points.
The rate of diffusion of substance i at mesh point k is
and thence
(11)
To determine completely the right-hand side of Eq. (11) is now necessary to find an expression for the activity coefficient, γ
i
, and the frictional coefficient, f
i
, contained in the expression of the diffusion coefficient. In fact, by substituting Eq. (2) into Eq. (7) we obtain the diffusion coefficient in terms of activity coefficients γ
i
:
(12)
where the frictional coefficient is assumed to be linearly dependent on the concentration of the solute like in sedimentation processes, i.e. in a mesh k, fi, k is
(13)
where k
f
is an empirical constant, whose value can be derived from the ratio R = k
f
/[η ]: Accordingly to the Mark-Houwink equation [27], [η ] = kMα is the intrinsic viscosity coefficient, α is related to the shape of the molecules of the solvent, and M is the molecular mass of the solute. If the molecules are spherical, the intrinsic viscosity is independent of the size of the molecules, so that α = 0. All globular proteins, regardless of their size, have essentially the same [η ]. If a protein is elongated, its molecules are more effective in increasing the viscosity and [η ] is larger. Values of 1.3 or higher are frequently obtained for molecules that exist in solution as extended chains. Long-chain molecules that are coiled in solution give intermediate values of is α, frequently in the range from 0.6 to 0.75 [27]. For globular macromolecules, R has a value in the range of 1.4 - 1.7, with lower values for more asymmetric particles [28].
Although Eq. (13) is a simplified linear model of the frictional forces, it works quite well in many case studies and can be easily extended to treat more complex frictional effects (see [25, 29]).
Let us focus now on calculation of the activity coefficients: a way to estimate the frictional coefficients will be presented in the next subsection. By using the subscript '1' to denote the solvent and '2' to denote the solute, it can be written that
(14)
where γ2 is the activity coefficient of the solute and c2 is the concentration of the solute. Differentiating with respect to c2 gives
(15)
The chemical potential of the solvent is related to the osmotic pressure (II) by
(16)
where V1 is the partial molar volume of the solvent and its standard chemical potential. Assuming V1 to be constant [30] and Differentiating μ1 with respect to c2 yield
(17)
Now, from the Gibbs-Duhem relation [31], the derivative of the chemical potential of the solute with respect to the solute concentration is
(18)
where M is molecular mass of the solute and is the partial molar volume of the solute divided by its molecular mass. The concentration dependence of osmotic pressure is usually written as
(19)
where B is the second virial coefficient, and thence the derivative of II with respect to the solute concentration is
(20)
Introducing Eq. (20) into Eq. (18) gives
(21)
From Eq. (15) and Eq. (21) it can be obtained that
so that
On the grounds that [32], solving the integral yields
(22)
The molecular mass M
i,k
of the species i in the mesh k can be expressed as the ratio between the mass, m
i,k
, of the species i in that mesh and the Avogadro number M
i,k
= m
i,k
/N
A
. If p
i
is the mass of a molecule of species i and c
i,k
· l is the number of molecules of species i in the mesh k, then the molecular mass of the solute of species i in the mesh k is given by
(23)
Substituting the expression in Eq. (22) gives, for the activity coefficient of the solute of species i in the mesh k (γ
i;k
), the following equation
(24)
Therefore, substituting Eq. (13) and Eq. (24) into Eq. (12), we obtain the following expression for a time- and space-dependent diffusion coefficient
(25)
We finally estimated in the following way the second virial coefficient B. The statistical mechanics definition of the second virial coefficient is as follows
(26)
where u(r), which is given in Eq. (27), is the interaction free energy between two molecules, r is the intermolecular center-center distance, k
B
is the Boltzman constant, and T the temperature. In this work, it is assumed that u(r) is the Lennard-Jones pair (12,6)-potential (Eq. 27), that captures the attractive nature of the Van der Waals interactions and the very short-range Born repulsion due to the overlap of the electron clouds:
(27)
By expanding the term exp into an infinite series, the Eq. (26) becomes
where T* ≡ 4/(k
B
T ) and thus
(28)
An estimate of B is given by truncating the infinite series of functions to j = 4, since simulation results not shown here prove that taking into account the additional terms, obtained for j > 4, does not significantly influence the simulation results [25].
Modelling the stochasticity
Both diffusion and reactions are modelled as reaction events whose dynamics is driven by the First Reaction Method of the Gillespie algorithm.
In particular, the diffusion events are modeled as first-order reactions. namely, the movement of a molecule A from box i to box j is represented by the reaction , where A
i
denotes the molecule A in the mesh i and A
j
denotes the molecule A in the mesh j. In this way, the reaction-diffusion system is modeled as a pure reaction system.
The space domain of the system is divided into a number N
s
of squared meshes of size l. The time evolution of the system is computed by the First reaction Method of Gillespie [24] that at each simulation step selects in each mesh the fastest reaction, compares the velocities of the N
s
selected reactions and finally executes the reaction that is by far the fastest. The fastest reaction is defined as the reaction whose waiting time is the smallest.
The time at which each event is expected to occur is a random variable extracted by an exponential distribution [24]. Let R
i
be the i-th reaction channel expressed as
where l
ij
is the stoichiometric coefficient of reactant S
p
(
i,j
), p(i, j) is the index that selects the species that participates in R
i
, L
i
is the number of reactants in R
i
, and r
i
is the rate constant. If the fundamental hypothesis of stochastic chemical kinetics [24] holds within a box, both diffusion and reaction events waiting times are distributed according to a negative exponential distribution, so that a typical time step has size
(29)
where R is the number of events. It is given by R = R
diff
+ R
react
, where R
diff
is the number of possible diffusion events and R
react
is the number of reaction events [33]. The diffusion and reaction propensities are given by the following expressions, respectively
(30)
(31)
where and are the number of chemical species that diffuse and the number of those the undergo to reactions, respectively. In general , since some species both diffuse and react. In Eq. (30), r
i
(diff)is the kinetic rate associated to the jumps between neighboring subvolumes, whereas in Eq. (31), r
i
(react)is the stochastic rate constants of the i-th reaction.
From Eq. (11), the rate coefficient of the first order reaction representing a diffusion event is recognized to be as follows
(32)