A Langevin dynamics approach for multi-layer mass transfer problems

Highlights

• Physically-based Langevin dynamics approach of diffusive mass transfer in a multi-layer medium.

• Conceptual and technical simple algorithm for simulating mass transfer through a sharp edge.

• Application to a drug-eluting stent model demonstrates the accuracy and usefulness of the method.

Abstract

We use Langevin dynamics simulations to study the mass diffusion problem across two adjacent porous layers of different transport properties. At the interface between the layers, we impose the Kedem–Katchalsky (KK) interfacial boundary condition that is well suited in a general situation. A detailed algorithm for the implementation of the KK interfacial condition in the Langevin dynamics framework is presented. As a case study, we consider a two-layer diffusion model of a drug-eluting stent. The simulation results are compared with those obtained from the solution of the corresponding continuum diffusion equation, and an excellent agreement is shown.

Introduction

Multi-layer diffusion problems arise in a number of applications of heat and mass transfer. Some industrial examples are moisture diffusion in woven fabric composites [1], hydrodynamics of stratified fluids and geological profiles [2], environmental phenomena such as transport of contaminants, chemicals and gases in layered porous media [3], and chamber-based gas fluxes [4]. Numerous applications concern the biomedical field and include, for example, transdermal drug delivery [5], drug-eluting stents [6] or brain tumor growth [7]. While here we focus on multi-layer diffusion, other related concepts such as anomalous diffusion, fractal kinetics and non-homogeneous layers, have been also studied within the context of drug release, see e.g., [8], [9], [10].

Often, the transported material is initially concentrated in one of the layers from which it propagates to the others by diffusion. The rate of transfer across the system in mainly determined by the diffusion coefficients in each layer. In many practical applications it is essential to regulate the mass flux between layers by suitable interface conditions. This can be accomplished, for instance, by placing a selective barrier between adjacent layers, which induces a chemical potential gradient at the boundary. Another mean for controlling the transfer rate are membranes, which are essentially very thin boundary layers with a small diffusion coefficient [11]. In addition to their role in slowing down the diffusion rate, membranes are also employed for specific functions, including separation/purification of gases, vapors, liquids, selection of ions, or other biological functions. Membranes are routinely used for medical care and individual protection, such as wound dressing, dialysis, tissue engineering, and controlled release of drugs. Membranes are also used for environmental cleaning and protection, such as water purification and air filtration. A better understanding of physical behavior of membranes as rate-controlling barriers can greatly improve the efficiency of separation and enhance their performance [12].

In this work, we consider simple models for mass transfer in multi-layered systems. We assume that the molecules are transported across the boundaries by passive diffusion only, i.e., no active transport process is performed to drive the random motion of molecules. Passive diffusion continues until enough molecules have passed from a region of higher to a region of lower concentration, to make the concentration uniform. When equilibrium is established, the flux of molecules vanishes: the molecules keep moving, but an equal number of them move into and out of both layers. Much work has been done from the analytical and computational point of view for treating multi-layer diffusion in continuum mechanics. An important aspect of layered systems is the matching conditions at the interfaces, where an interface is the common boundary between two layers. Analytical solutions to such problems are highly valuable as they provide a great level of insight into the diffusive dynamics, and can be used to benchmark numerical solutions [13]. Various methods are available for the analysis and the solution of such problems [14], [15]: The orthogonal expansion technique and the Green’s function approach [16], [17], [18], [19], the adjoint solution technique [20], the Laplace transform method [14], [15], [21], [22], [23], and finite integral transforms [24], [25], [26]. Integral transform techniques applied to heat transfer problems was discussed in great detail in the book by Özişik [20], where several different transformations are given depending on the situation. However, there are severe numerical instabilities and computational drawbacks that arise when the number of layers increases [22]. Other papers demonstrate the complexity of solving diffusion problems with a large number of layers, either using eigenfunction expansion for somewhat different boundary conditions [27], or based on the Green function approach with biological applications [28]. Computational complexity of finite difference schemes is widely discussed [29].

Recently, a new computational method for studying diffusion problems in multi-layer systems has been proposed [30], [31]. The method is based on the well-established notion that Brownian dynamics of particles can be also described by the Langevin’s equation (LE) [32]. Therefore, the particle’s probability distribution function (or, equivalently, the material concentration) can be computed from an ensemble of statistically-independent single particle trajectories generated by numerical integration of the corresponding LE. Integrating LE within each layer is pretty straightforward, and there are a number of algorithms (Langevin “thermostats”) that are widely used for molecular dynamics simulations at constant temperature [33], [34], [35]. The key problem is how to perform the integration during time-steps where the particle moves between layers, in a manner ensuring that the imposed interlayer conditions are satisfied. In Ref. [31], a set of algorithms for handling the dynamics across sharp interfaces has been introduced. Here we present an algorithm that combines many types of interfaces (a sudden change in diffusivity, a semi-permeable membrane, an imperfect contact), with the advantage of treating all these cases with a unified physical-based method. The new algorithm is applied for studying a two-layer model of drug release from a drug eluting stent into the artery. Excellent agreement is found between the LE computational results and the semi-analytical solution.