Ratcheted diffusion transport through crowded nanochannels

Abstract

The problem of transport through nanochannels is one of the major questions in cell biology, with a wide range of applications. In this paper we discuss the process of spontaneous translocation of molecules (Brownian particles) by ratcheted diffusion: a problem relevant for protein translocation along bacterial flagella or injectosome complex, or DNA translocation by bacteriophages. We use molecular dynamics simulations and statistical theory to identify two regimes of transport: at low rate of particle injection into the channel the process is controlled by the individual diffusion towards the open end (the first passage problem), while at a higher rate of injection the crowded regime sets in. In this regime the particle density in the channel reaches a constant saturation level and the resistance force increases substantially, due to the osmotic pressure build-up. To achieve a steady-state transport, the apparatus that injects new particles into a crowded channel has to operate with an increasing power consumption, proportional to the length of the channel and the required rate of transport. The analysis of resistance force and accordingly – the power required to inject the particles into a crowded channel to overcome its clogging, is also relevant for many microfluidics applications.

Introduction

Ratcheted diffusion can be described as directionally biased random thermal motion. Physically, this process is often described as a transport phenomenon in spatially periodic systems, a typical example of which is a “saw-tooth” potential that is periodic in space and has a broken forward-backward symmetry, however, ratchet potentials exist in various forms in terms of symmetry and periodic behavior¹. The function of biologically relevant “motor” enzymes, or molecular motors, is based on non-equilibrium ratcheted diffusion, where the energy is continuously recycled in the system – typically by the ATP hydrolysis. Well studied classes of these molecular motors are myosins, kinesins and dyneins¹. The physical mechanism producing the powered linear motion along the filament with asymmetric (saw-tooth) potential, energized by the ATP hydrolysis, is well understood after the classical work of Prost et al.^2,3. Myosins move along actin filaments and these proteins are responsible for muscle contraction⁴. Kinesin and dynein motor proteins move along microtubules, either towards the elongating/polymerising + end (which has a high concentration of so called β-subunits) or the − end (which has a high concentration of α-subunits). Dyneins have various additional functions which include facilitating the movement of cilia and flagella⁴.

Our interest is different: we examine the equilibrium situation when the directed transport is achieved spontaneously, energized by the Brownian thermal motion and directed by restricting the diffusion in the “wrong” direction. Examples of this equilibrium diffusive transport is the translocation of cargo across membranes through nanopores or along nanochannels⁵. Essentially, the directed Brownian motion is achieved by simply not allowing the cargo to travel back to the side of the channel where the process of translocation started. To achieve this, the cargo has to block the channel width such that no overtaking is possible in a single-file (one-dimensional) motion⁶. In this case the potential profile along the channel is irrelevant, as the equilibrium detailed balance would not produce any directed motion without an external force (in spite of any symmetry breaking). A variety of chemical asymmetries may exist inside the channel to affect the movement of the cargo⁷. However, the effect of the channel potential is merely to modify the effective diffusion constant, which is again a well-studied problem⁸. Biological examples of transport assumed to rely on spontaneous ratcheted diffusion are found in various biological systems, such as:

• 1 flagella – hair-like protrusions from the cell body of specific pro- and eukaryotic cells⁹ grow by transporting flagellin protein subunits from the body of bacteria through a 2 nm wide channel over distances of 10–20 μm, at a constant rate¹⁰;

• 2 the injectosome complex (type III secretion system) – a needle-like protein appendage that is found in gram-negative bacteria, the function of which is to secrete proteins into eukaryotic host cells to facilitate the process of infection¹¹;

• 3 the general secretory (Sec) pathway, in which the translocase protein complex machinery provides a translocation pathway for proteins from the cytosol across the cytoplasmic membrane in bacteria¹²;

• 4 bacteriophage (bacterial virus) DNA injection into host bacterial cells through the tail of the phage – a good example of this process is the T4 phage mechanism of injection¹³.

Independently, ratcheted transport in Brownian systems is of great practical relevance for applications in lab-on-chip microfluidic devices¹⁴. Typically, if crowding conditions of drops/bubbles/particles are achieved at the inlet of the microfluidic channel, the resistance to the transport develops, which eventually leads to microchannel clogging¹⁵. The latter is a practically important phenomenon which limits the performance of microfluidic devices and of industrial micromixers^16,17.

One particular biological system that we shall frequently refer to is the bacterial flagellum, which is a tube with a 2 nm wide inner channel and the walls formed by crystallized flagellin proteins. In order for the flagellum to grow, the protein subunits have to be translocated to the distal (growing) end of the flagellum through this channel, which is too narrow to accommodate folded proteins. Hence, subunits are transported in their unfolded state, after they are inserted into the channel by the export apparatus¹⁸. Flagellar export apparatus is a multicomponent complex, energized by ATP hydrolysis and achieving the successive unfolding of flagellin chains and insertion of them into the channel^19,20,21, thus providing a bias of the subsequent 1D Brownian motion, which motivates the work presented in this paper.

Very recently there has been an attempt on studying the diffusion of flagellin subunits along this channel²², where many of the details and parameters of this biological system are carefully presented. The shortcoming of this basic simulation is that it does not take into account interactions (essentially – repulsion on contact) between travelling particles, which then fails to describe the crowding regime. Here we show molecular-dynamics simulation results for translocation of Brownian particles in nanochannels under conditions at the inlet, which mimic crowded cellular environments. The effect of molecular crowding at the channel inlet is twofold: on the one hand it acts as a net drift speeding up the translocation near the open outlet, but on the other hand it also causes an increasing resistance force, slowing down the transport near the channel inlet. The range of channel length where crowding contributes a net resistance force increases with the total channel length. Remarkably, this effective resistance force is an emergent phenomenon which in 1D depends entirely on the cooperative multi-body dynamics of particles colliding in the channel. Collective diffusion combined with a third-order virial expansion provides a simple theoretical description of the density profile inside the channel in agreement with the simulation data. The transition from classical free diffusion to the crowding regime is controlled by the competition between the particle injection rate and the mean first passage time of free diffusion. This fact is crucial in order to assess the onset of crowding in both biological and engineering systems.

Results

Free single-file diffusion

The ratcheting, or rectification of the Brownian motion along the narrow channel is achieved by continuously inserting new particles from one end, thus restricting the backward motion – while the channel remains open at the distal end. The diffusion process as described in this paper raises three main issues: 1) how high is the resistance force, if the particles are inserted at a sufficiently high rate? 2) How do the two characteristic times, the mean first passage time τ_diff and the rate of insertion 1/T_in compete? 3) What is the crossover between the stochastic diffusion dominated regime (when ) and the dynamically restricted, resistance-dominated regime in a crowded channel?

In order to introduce the problem, we first look at the motion of a single particle inside a channel blocked at one end, in order to verify the properties of generalized diffusion, fig. 1. Since the particle essentially disappears as soon as it reaches the open end, a distance L away, this is a problem of first passage with a reflecting boundary condition at x = 0 and an absorbing condition at x = L. The probability distribution of first passage times is a classical result obtained by the method of images due to Lord Kelvin²³; it takes the form (see “Methods” for detail):

This probability distribution is used to fit the histogram of simulation results for the frequency of diffusion time of a single particle along the channel, obtained by the Brownian dynamics simulation, fig. 1. In particular, the free diffusion constant in our simulation is thus measured to be D = 0.0065 in the reduced simulation units of [σ²/τ]. An example analysis in “Methods” shows that this corresponds to a value D = 6.5 · 10⁻¹² m²/s for 1-micron size polystyrene colloids in water.

Figure 1

A histogram representing the distribution profile of individual diffusion times, τ_diff, along the channel of length L = 60σ (σ being the size of particle).

The time is measures in [kts]: 1000 of simulation time steps. The solid line represents the fit to the probability distribution discussed in the text, eq.(1), with only the normalization factor and the diffusion constant D as fitting parameters.

Unlike earlier work^1,8, we do not consider designed periodic potentials in the channel for simplicity and due to the fact that this is not necessary in order to find D, the diffusion constant. The important result, reproduced frequently in the literature, is that the effective diffusion constant for traveling in a channel with a potential V(x) along its length is related to the bare diffusion constant (with no potential) as²⁴

The product of two integrals in the denominator is always greater than unity (this is demonstrated by applying the Cauchy-Schwartz inequality²⁴, see also²⁵) and therefore, for any potential V(x) on the channel walls the diffusion is hindered. For our purposes of studying the effects of crowding in the channel, we may simply take a given fixed value for the effective diffusion constant, D, obtained by the fitting of single-particle diffusion data as in fig. 1 and proceed to examine the case when many particles restrict their motion by crowding.

Crossover to crowded regime

Having determined the average time of individual particle diffusion to the far end of the channel, τ_diff = ∫t · f(t)dt = L²/D, we can now simulate the process with particles injected at x = 0 with increasingly short time intervals T_in (i.e. at an increasing rate, or the particle flux J_in = 1/T_in). As explained in “Methods” below, our simulation algorithm only allows the insertion if a sufficient room (of the particle size σ) is available at x = 0; if the space is occupied, the insertion is missed and is attempted after the next T_in interval. At low rate of insertion this omission never happens, but at this becomes a limiting factor. Figure 2 represents the results in the form of the cumulative number of particles having passed through the channel: N(t) is essentially the sequential number of a particle leaving the channel at a time t. From the plot it is clear that at sufficiently long times of simulation the progression of particles is linear on average. One may interpret this as a constant rate of delivery of cargo along the channel, however, we shall discuss below the constraints on the required power and the force exerted by the insertion apparatus.

Figure 2

A graph of a cumulative number of particles having passed through the channel of L = 30σ against simulation time-step (measured in [kts]: 1000's of time steps, as in fig. 1).

The labeled values of T_in, also in [kts], represents the set interval of an attempted insertion at x = 0. The large scale of this cumulative plot makes obvious the linearity of N(t), but masks what happens at shorter times – when the insertion location may be occupied by other particles; fig. 3 below illustrates the effect of crowding restriction.

Figure 3 plots the average slope of the data in fig. 2 on the double-logarithmic scale to allow examination of the two key findings. We choose to plot the inverse flux, dt/dN(t) against T_in, to have a more clear sense of the events. The plot also shows the results for three different channel lengths (measured in the units of particle size σ), as labeled on the figure. In essence, every curve in fig. 2 gives a single point in this graph, for a given L and T_in. Two regimes are evident: when particles are inserted at a sufficiently low rate, the insertion rate is simply equal to the constant flux J₀ through the channel, with a small correcting factor reflecting the channel length. However, when the attempted rate of insertion becomes high, the rate of diffusive transport through the channel saturates at a constant maximum value that depends on the channel length and is determined by the increasing resistance for the motion of the single file of particles packed in the channel.

Figure 3

The average slope of the time-step versus the cumulative number of particles as a function of insertion time, T_in for three channel lengths.

For high T_in, Δt/ΔN = T_in. The slope remains constant as a function of T_in initially due to the fact that insertion rate is too high and particles cannot be inserted more often than T_in ≈ 1000 ts, the crossover between the constant insertion rate and the linearly increasing rate at higher T_in.

Since we are interested in the problem of crowding, the most valuable scenario arises when the insertion rate is high and so particles frequently collide in the channel. In order to understand how crowding inside the channel affects the interactions between particles, we considered two scenarios – an open channel, as described earlier, in which the particles essentially disappear at the far end (i.e. not able to return back into the channel once they reach x = L) – and a closed channel with another reflecting wall at x = L. The latter channel becomes filled with particles at a specific rate. Because the length of the channel emerges as an important parameter, we have again compared channels of three lengths: L = 30σ, 60σ and 90σ, fig. 4. As the time increases, the particles initially fill the channel at a constant rate – but then their number in channel saturates at a constant value: this is a lower and a highly fluctuating value if the channel is open and a higher and almost constant value if the channel is closed and the particles pack until there is no more insertion possible. But even the latter value is much less than the channel length: 79 in the channel of 90-particles long, 55 in 60 and 28 in 30. The difference is due to the effective resistance to insertion of new particles when the density increases sufficiently; note how this difference grows with the length of the channel, which reflects the growing resistance force for new particle insertion.

Figure 4

The number of particles in open and closed channels of different lengths: L = 30σ, 60σ and 90σ for the fastest attempted insertion rate (T_in = 0.001 kts); the curves are labeled accordingly on the plot.

The fluctuations in open channels are caused by the diffusive flux of particles exiting the channel. Due to a higher overall number of particles in longer channels, the osmotic pressure inside is higher, which explains the larger gap between L and the mean plateau values for longer channels.

The growth of number of particles in the channel as a function of time was fitted with the simple exponential function: n(t) = n_max(1 − exp[−t/τ_rel]), where n_max is the maximum number of particles in a channel at a given time (fixed for a channel of given length): the mean plateau value in fig. 4 gives n_max = 54 for L = 90σ, 34 for 60σ and 17 for 30σ. The gap between the dense-packing and the plateau value in this case (of an open channel) grows approximately linearly with the channel length, in contrast with the closed channel where this gap increases approximately as the square of length. Remarkably, the characteristic relaxation time τ_rel after which the dense-packing plateau is achieved was the same in all cases, τ_rel ≈ 40 kts.

Resistance to transport

Here we show that the effect of crowding inside the channel results in an effective potential within the channel when the density of particles ρ is high. Let us consider a steady state, at t > τ_rel, when the flux is constant along the channel, J_in = J_out. The simple one-dimensional diffusion with the boundary conditions (ρ₀ at x = 0 and ρ = 0 at r = L) gives the time-evolving profiles which eventually would reach the steady state of linear ρ(x) = [J_eq/D](L − x) and constant flux J_eq = −D∇ρ. This steady state is only possible if the density at the entrance ρ(x = 0) = [J_eq/D]L remains small, otherwise crowding becomes relevant in the problem and the simple diffusion solutions no longer valid.

In order to examine crowding regime in more detail, fig. 5 shows the profiles of particle density within the channel at steady-state (in other words, when the channel is filled at t > τ_rel). The density ρ(x) here is defined as the number of particles per unit length, measuring the probability density of finding a particle at a position x in the channel. The deviation from the linear density profile is prominent. Generalized diffusion can be used to describe transport at arbitrary particle concentrations^26,27. The relevant equation reads:

where D_coll = D∂(βΠ)/∂ρ, with β = 1/k_BT, is the collective diffusion coefficient, as defined in nonequilibrium thermodynamics²⁸. It accounts for deviations from free diffusion through the osmotic pressure Π(ρ), which for a non-ideal fluid can be written as a virial expansion in the density βΠ = ρ + B₂ρ² + B₃ρ³, with B₂ and B₃ the second and third virial coefficients, accounting respectively for two-body and three-body interactions between particles in the channel. In bulk 3D systems, higher order coefficients are usually needed in the expansion to properly describe crowding effects at high density²⁷, however, we shall see below that in our case of single-file 1D transport terms up to third order are sufficient to capture the phenomenology observed by simulations. This is simply because one particle in front and one behind, determine the progress of any given particle. The steady state solution of eq.(3) is easy to obtain by the chain rule:

The osmotic pressure has to be a linear function of coordinate along the channel! When the density is low and Π = k_BTρ, this reflects the classical steady-state diffusion with the linear density variation ρ(x). At high density we need to invert the dependence Π(ρ). Within the third-order virial approximation we obtain

where the virial coefficients B₂ and B₃ are both positive for repulsive particles used in our simulation. The coefficients b and c used as fitting parameters to successfully reproduce the measured density profiles in fig. 5 are explicitly expressed here and we see that b < 0 always, while the sign of c depends on the relation between B₂ and B₃. For our repulsive Lennard-Jones potential the fitting suggests that .

Figure 5

The steady-state local density of particles along the channel of varying lengths (labeled on plot) when it is filled at the maximum rate T_in = 0.001 [kts].

Error bars represent 1 standard deviation based on 80 snapshots at different times during the equilibrium stage of simulation. The solid lines represent theoretical fitting curves obtained with ρ = a(L − x) + b(L − x)² + c(L − x)³ (see eq.(5)), where a = (0.84/L − bL − cL²) and b = −0.0014, −0.0007, −0.0004, c = 1.8 · 10⁻⁵, 6.0 · 10⁻⁶, 2.2 · 10⁻⁶ for L = 30,60,90σ, respectively.

The eq.(4) for the osmotic pressure illustrates the first barrier to transport into such a channel. The export apparatus that injects new particles into the channel, at x = 0, has to exert a force F = −(1/ρ)[dΠ/dx] to overcome a pressure gradient at the channel entrance: Π(x = 0) over a distance comparable to the particle size σ. This “injection force” is F_in ≈ k_BT[J_eq/D]L, linearly increasing with the channel length if the particle flux is to be maintained constant, as in¹⁰. Accordingly, the power required to maintain the constant flux will also be proportional to the channel length: .

The effect of crowding on the Brownian resistance force along the channel is described in greater detail in the Methods section. There we show that when the two aditional virial-coefficient terms make the steady-flux density distribution nonlinear, eq.(5), the local force acting on the particles along the channel acquires an additional contribution due to the collective Brownian effects:

This equation predicts that the (negative) drag force is linear in L for suffiently long channels and hence proportional to the total number of particles in a long channel.

Discussion

We have studied ratcheted diffusion transport in nanochannels where translocation is achieved by means of crowding at the channel inlet and depletion at the outlet. This problem is paradigmatic for biomolecule translocation in cellular pores and flagella⁹, as well as for microfluidics where crowding leads to microchannel clogging^15,16. In a nanochannel, the cargo (protein subunits in a flagellar channel or colloids in a microfluidic channel) has a size comparable to the channel width: the particles move in a single file towards the distal open end, after being injected by an appropriate export machinery from the proximal end. We find many interesting aspects of such motion when the density of particles in the channel becomes high, i.e. their mutual collisions start having a significant effect on the diffusion (crowded regime). Of a particular importance is the criterion of when the crowding regime sets in, which is determined by comparing the injection period T_in = 1/I_in and the time of diffusion along the channel, τ_diff = L²/D. When the insertion rate is high enough, it takes a certain time to get into the crowded regime; this saturation (clogging) time appears to be independent of the channel length, because it is against the resistance force to the column motion that the clogging occurs.

But the main message and conclusion are expressed at the end of the last section: in order to propel the long column of particles along the channel, maintaining the constant (steady-state) current, the inserting apparatus has to exert an increasing force and expend an increasing power – both approximately linearly increasing with the channel length. For long enough channels, the motion will simply stop (channel clogging) until the density slowly drops after gradual ejection of particles from the open end and redistribution of remaining ones along the length. For biological systems such as the flagellin transport, the ATP-powered export apparatus would not be able to increase its power output with the growth of the flagellar length and instead would have to slow down the rate of insertion to compensate for the growth of resistance. Empirical evidence so far points towards a constant growth rate of flagella and further investigations are required in the future to clarify the actual mechanism.

In crowded environments, non-specific interactions between molecules lead to deviations from free diffusion^29,30. Here we have shown that repulsive interactions between particles lead to an effective drag force, slowing down diffusive transport, which is an emergent cooperative phenomenon. This is an important effect which has to be taken into account in the analysis of many both natural and engineered systems: from biomolecule translocation in membrane pores to the design of microfluidic lab-on-chip devices.

Methods

We approach the problem of ratcheted diffusion by Brownian dynamics simulations. Our setup consists of a cylindrical channel into which particles are added at a certain rate. The force field of the model is rather simple – there is only hard core repulsion between particles as well as between particles and the wall of the cylindrical channel, which is represented by the repulsive part of the Lennard-Jones potential. The channel is blocked from one side, which only allows particles to translocate to the other side of the channel. As soon as particles reach the end of the channel they disappear from the system. The dynamics of the process of translocation as well as particle's positions and forces acting and each particle were recorded at every time-step. The number of injected - as well as ejected particles was also monitored. Based on this information, the diffusion profile inside the channel can be defined as soon as the injection/ejection equilibrium state is reached in the channel, meaning that at this point, the channel is maximally occupied and the rate of injection is equal to the rate of ejection. New particles are only added when the position at the beginning of the channel is not occupied by other particles and therefore even if the particles are set to be added every time-step, the addition will not happen if another particle is found less than σ distance away from the addition site. Additionally, we studied how the system behavior changes if the channel is closed, as opposed to deleting particles when they reach the end of the channel.

In order to be able to convert the results of simulations, Lennard-Jones (reduced) units were used. In molecular dynamics algorithms, physical quantities are expressed as dimensionless units, which results in all physical quantities being around unity, hence reducing the problems of working with numbers that are very small or very large, making errors more visible and also simplifying the output produced and increasing computational efficiency as the equations of motion become simplified due to the absorption of model-defining parameters into reduced units. Reduced units also allow one to work on different problems with a single model. The following reduced units are used for a Lennard-Jones system: length - σ, energy - ε and mass m. From these, one can express units of all other physical quantities as follows: time (τ) as , temperature k_BT in units of ε and pressure P in units of ε/σ³. The simulation time step was taken as ts = 0.01τ.

For instance, a ‘typical’ colloid particle of size σ = 1 μm and mass density 1000 kg/m³ would have a mass m ≈ 4.2 · 10⁻¹⁵ kg (staying strictly in SI units here). Taking k_BT = ε ≈ 4.2 · 10⁻²¹ J (for room temperature), we obtain the time scale τ = 10⁻³ s. This means that the diffusion constant that we obtained from the fit in fig. 1, D ≈ 0.0065 in units [σ²/τ] would have the value D ≈ 6.5 · 10⁻¹² m²/s, which is almost exactly the magnitude of a widely measured diffusion constant for micron-size sterically stabilized polystyrene particles in water.

First passage time

Here we summarize the standard facts about the first passage problem. This is the classical problem addressing the question: if a diffusing particle starts at, e.g. x = 0 – on average how long would it take to reach a given point, x = L, for the first time (hence the name of “first passage”). Obviously in the stochastic problem one can only have the answer for the average time and so the aim is to find the probability distribution of such times (and then evaluate the average if required). In our case we remain in the 1-dimensional limit and the generic diffusion equation for the concentration, or the probability to find a particle at (x,t), has the form of continuity relation:

where J(x,t) is the generalized diffusion flux. For the free diffusion, the Green's function of this equation with is the Gaussian

The “first passage” condition is implemented by applying the absorbing boundary condition, that is, once the particle first reaches x = L, it disappears from the analysis: p(L,t) = 0. Such a condition is traditionally treated by the method of images, that is, imposing two opposite values – one starting from x = 0, the other from x = 2L:

Figure 6 plots several curves of this probability distribution, for D = 1 and the increasing time (the particles start from x = 0 and x = 2L at t = 0), to show how the method of images helps to maintain the constraint of p(L) = 0 at all times.

Figure 6

The probability distribution Eq.(9) for particles with the absorbing boundary enforced at x = L by the method of images.

The survival probability for the particle to remain anywhere between 0 and L, having started at x = 0 and t = 0 is obtained by integration: . The result of this integration is the difference of two error functions:

Since x = L is an absorbing boundary, the domain of 0 − L is independent from the domain L − 2L. Given the definition of the survival probability, the fraction −∂Q(t)/∂t is absorbed between t and t + dt. This means that f(t) = −∂Q(t)/∂t is actually the probability density of the time t that takes the particle to reach x = L for the first time, given in the main text as eq.(1). Figure 7 plots a collection of these distributions for different values of the diffusion constant and fixed L = 1. In both plots we did not normalize the probability density functions or their parameters, just assuming all values are non-dimensional for qualitative illustration. To find the average value of the first passage time we need to integrate, obtaining the very much expected result:

Figure 7

The probability distributions eq.(1) of the first passage times in the channel of 0 − L, for several values of the diffusion constant D.

Forces in crowded channel

The effect of crowding on the ratcheted transport can be understood via the Brownian force acting on a particle in the system is given in general terms as . Since we work in 1-dimensional geometry, with particles unable to overtake each other, only the horizontal projections of all forces matter. We can take as positive all forces parallel to the x axis along the translocation direction (down the density gradient) and with the minus sign we take the resistance forces which oppose the translocation, i.e. drag forces. Using the chain rule,

and employing the virial expansion for Π(ρ) and expression for ρ(x) given by the eq.(5) of the main text, after minor manipulation we obtain the additional effective force determined by the collective effects in the crowded regime:

where the virial coefficients B₂ and B₃ were inroduced in the main text. The last term can be negative, i.e. contribute a net drag force or resistance to diffusive motion. It is easy to carry out the explicit calculation of the total drag force per particle due to crowding in the channel, by taking the integral over the channel length, with the result given in eq.(6).

It is also interesting to observe that the effect of crowding on diffusive transport can be mapped onto a classical diffusion equation in the effective force-field given by F_crowd, in the spirit of an approximation introduced in²⁷,

This is a classical (steady-state) diffusion equation in the force-field F_crowd which takes care of the additional ρ-dependent crowding effect. Hence, the diffusive transport in a crowded environment can be effectively described by means of classical diffusion in the additional force-field due to crowding. The latter can be calculated analytically using (12) whenever the inter-particle interaction and the corresponding virial coefficients are known. The solution to this equation is given by the expression for ρ that we presented as eq.(5).