A Langevin dynamics approach for multi-layer mass transfer problems
AA Langevin dynamics approachfor multi-layer mass transfer problems
Oded Farago and Giuseppe Pontrelli ∗ July 28, 2020
Abstract
We use Langevin dynamics simulations to study the mass diffusion problem across twoadjacent porous layers of different transport property. At the interface between the layers,we impose the Kedem-Katchalsky (KK) interfacial boundary condition that is well suitedin a general situation. A detailed algorithm for the implementation of the KK interfacialcondition in the Langevin dynamics framework is presented. As a case study, we consider atwo-layer diffusion model of a drug-eluting stent. The simulation results are compared withthose obtained from the solution of the corresponding continuum diffusion equation, and anexcellent agreement is shown.
Keywords : composite materials, interface conditions, diffusion equations, mass flux, Langevindynamics
Multi-layer diffusion problems arise in a number of applications of heat and mass transfer. Someindustrial examples are moisture diffusion in woven fabric composites [1], hydrodynamics of strat-ified fluids and geological profiles [2], environmental phenomena such as transport of contami-nants, chemicals and gases in layered porous media [3], and chamber-based gas fluxes [4]. Numer-ous 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 dif-fusion, other related concepts such as anomalous diffusion, fractal kinetics and non-homogenouslayers, have been also studied within the context of drug release, see e.g., [8, 9, 10]. ∗ Corresponding author. Email: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] J u l ften, the transported material is initially concentrated in one of the layers from which it prop-agates to the others by diffusion. The rate of transfer across the system in mainly determined bythe diffusion coefficients in each layer. In many practical applications it is essential to regulate themass 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 gradientat the boundary. Another mean for controlling the transfer rate are membranes which are essen-tially very thin boundary layers with a small diffusion coefficient [11]. In addition to their role inslowing down the diffusion rate, membranes are also employed for specific functions, includingseparation/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 dress-ing, dialysis, tissue engineering, and controlled release of drugs. Membranes are also used forenvironmental cleaning and protection, such as water purification and air filtration. A better un-derstanding of physical behaviour of membranes as rate-controlling barriers can greatly improvethe efficiency of separation and enhance their performance [12].In this work, we consider simple models for mass transfer in multi-layered systems. We as-sume that the molecules are transported across the boundaries by passive diffusion only, i.e., noactive transport process is performed to drive the random motion of molecules. Passive diffusioncontinues until enough molecules have passed from a region of higher to a region of lower concen-tration, to make the concentration uniform. When equilibrium is established, the flux of moleculesvanish: 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 matchingconditions at the interfaces, where an interface is the common boundary between two layers. An-alytical solutions to such problems are highly valuable as they provide a great level of insight intothe diffusive dynamics and can be used to benchmark numerical solutions [13]. Various methodsare available for the analysis and the solution of such problems [14, 15]: The orthogonal expan-sion 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 reported in great detail in thebook by ¨Ozis¸ik [20], where several different transformations are given depending on the situation.However, there are severe numerical instabilities and computational drawbacks that arise whenthe number of layers increases [22]. Other papers demonstrate the complexity of solving diffu-sion problems with a large number of layers, either using eigenfunction expansion for somewhatdifferent boundary conditions [27], or based on the Green function approach with biological ap-plications [28]. Computational complexity of finite difference schemes is widely discussed [29].Recently, a new computational method for studying diffusion problems in multi-layer systemshas been proposed [30, 31]. The method is based on the well-established notion that Browniandynamics 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 becomputed from an ensemble of statistically-independent single particle trajectories generated bynumerical integration of the corresponding LE. Integrating LE within each layer is pretty straight-forward, and there are a number of algorithms (Langevin “thermostats”) that are widely used formolecular dynamics simulations at constant temperature [33, 34, 35]. The key problem is how to2erform the integration during time-steps where the particle moves between layers, in a mannerensuring that the imposed interlayer conditions are satisfied. In Ref. [31], a set of algorithms forhandling the dynamics across sharp interfaces has been introduced. Here we present an algo-rithm that combines many types of interfaces (a sudden change in diffusivity, a semi-permeablemembrane, and an imperfect contact), with the advantage of treating all these cases with a unifiedphysical-based method. The new algorithm is applied for studying a two-layer model of drugrelease from a drug eluting stent into the artery. Excellent agreement is found between the LEcomputational results and the semi-analytical solution. Let us consider a composite medium consisting of a number of layered slabs. A slab is definedhere as a plate that is homogeneous and isotropic, having a finite thickness, but extends to infinityin the other two dimensions. In a typical diffusion problem driven by concentration gradient,most of the mass dynamics occurs along the direction normal to the layers. We, therefore, restrictour study to a simplified one-dimensional model across a multi-layer system. The concentrationof material in each region, c i ( x, t ) ( i = 1 , . . . , n ), is governed by the time-dependent diffusionequation ∂c i ∂t = D i ∂ c i ∂x , (2.1)where D i is the diffusion coefficient in the i -th region. The concentrations in the adjacent regions i and i + 1 must be matched at the boundary between them, which is located at x = L i . Twointerfacial boundary conditions (IBCs) must be specified at each interface. If mass is conserved(no source or sink) at the interface, then the concentration flux must be continuous J i = − D i ∂c i ∂x = − D i +1 ∂c i +1 ∂x = J i +1 at x = L i , t > . (2.2)The other IBC to be specified at x = L i depends on the nature of the interface. The transportof material can be completely blocked by placing a perfectly reflecting ( J i = 0 ) or perfectlyabsorbing ( c i = 0 ) barriers. Typically, however, we are interested at intermediate situations wherethe mass flux is not completely blocked, but only hindered by interfaces whose aim is to controlthe rate of mass transfer across the layers. Here, we consider Kedem-Katchalsky (KK) IBC thatreads [36, 37] J i = P i ( c i − σ i c i +1 ) , at x = L i , t > , (2.3)where P i and σ i are, respectively, the permeability and partition coefficients of the KK condition.We focus on the KK IBC (2.3) because it represents the most general case of an interface whereboth a discontinuity in the chemical potential and a semi-permeable membrane are present, inaddition to a possible discontinuity in the diffusion coefficient. The case without a membranecorresponds to the limit P i → ∞ , when the KK IBC must be replaced with c i = σ i c i +1 , at x = L i , t > (2.4)3r, otherwise, the flux diverges at the interface. Eq. (2.4) describes the interfacial condition atan imperfect contact boundary with partition coefficient σ i arising from the discontinuity in thechemical potential of the transported molecules in the adjacent layers [31]. In the special case ofEq. (2.4) when σ i = 0 (or, σ i → ∞ ), we have c i = 0 (or, c i +1 = 0 ), which describes a perfectlyabsorbing boundary. A subcase of (2.4) is σ i = 1 (a perfect contact), when the concentrationexhibits no discontinuity for P i → ∞ . However, when P i is finite in eqn (2.3), we expect aconcentration jump even for σ i = 1 , as the KK IBC reduces to J i = P i ( c i − c i +1 ) , at x = L i , t > , (2.5)which is the IBC describing the effect of a thin semi-permeable membrane with permeability P i , but without a chemical potential jump. Finally, when P i = 0 , we recover the condition at aperfectly reflecting boundary, J i = 0 . The method presented in ref. [31] is based on the description of the overdamped Brownian motionof particles via the underdamped LE m dvdt = − α ( x ) v + β ( t ) + f ( x ) , (3.1)where m and v = dx/dt denote, respectively, the mass and velocity of the diffusing particle.This is Newton equation of motion under the action of a “deterministic” force f ( x ) . The im-pact of the random collisions between the Brownian particle and the molecules of the embeddingmedium is introduced by two additional forces - (i) a friction force, − α ( x ) v , and (ii) stochasticGaussian thermal noise, β ( t ) , with zero mean, (cid:104) β ( t ) (cid:105) = 0 , and delta-function auto-correlation, (cid:104) β ( t ) β ( t (cid:48) ) (cid:105) = 2 k B T α ( x ( t )) δ ( t − t (cid:48) ) , where T is the temperature and k B is Boltzmann’s con-stant [38]. The friction coefficient, α , in LE and the diffusion coefficients, D , in the correspondingdiffusion equation, satisfy the Einstein’s relation [32, 39]: α ( x ) D ( x ) = k B T. (3.2)In the Langevin dynamics approach to multi-layer diffusion, the concentration profile, c ( x, t ) ,is computed by generating an ensemble of statistically-independent particle trajectories of duration t , from which a fine-grained histogram can be constructed. We define c ( x, t ) such that, at t = 0 ,the total density is normalized to unity and essentially represents the initial probability distributionfunction of the particles (cid:90) ∞−∞ c ( x, dx = 1 . (3.3)The trajectories are calculated by numerically integrating Eq. (3.1). To allow for simulations ofLangevin dynamics in multi-layer systems, algorithms were derived in [31] for handling the tran-sition in presence of (i) layers with different diffusion coefficients, (ii) a semi-permeable mem-brane, and (iii) a step-function chemical potential. Here, we integrate them into a single unified4 igure 1: The typical two-layer one-dimensional system. A continuous flux [see Eq. (3.5)] is imposed at the interface x = 0 (dashed line), together with Kedem-Khatchalsky (KK) condition [see Eq. (3.6)]. algorithm for crossing a KK IBC [Eq. (2.3)] with continuous flux [IBC Eq. (2.2)]. We will notrepeat the discussion on the physical basis underlying the method, but rather present a practical recipe describing how to implement the algorithm. To this purpose, we consider the two-layersystem shown in fig. 1, with a step diffusion function D ( x ) = D x < D x > . (3.4)The continuity of flux J applies at the interface − D ∂c ∂x = − D ∂c ∂x at x = 0 , t > , (3.5)together with the KK IBC J = P ( c − σc ) , at x = 0 , t > . (3.6) The initial position of the particle is drawn from the probability distribution c ( x, , Eq. (3.3), andthe initial velocity from the Maxwell-Boltzmann distribution ρ MB ( v ) = (cid:114) m πk B T exp (cid:18) − mv k B T (cid:19) . (3.7)The trajectory x ( t ) is then computed by performing discrete-time integration of LE (3.1). For thispurpose, we use the algorithm of Grønbech-Jensen and Farago (GJF) [35] x n +1 = x n + b (cid:20) dtv n + dt m f n + dt m β n +1 (cid:21) (3.8) v n +1 = a v n + dt m (cid:0) a f n + f n +1 (cid:1) + bm β n +1 , (3.9)5o advance the coordinate x n = x ( t n ) and velocity v n = v ( t n ) by one time step from t n = n dt to t n +1 = t n + dt . In the above GJF equations (3.8)-(3.9), f n = f ( x n ) , β n is a Gaussian randomnumber satisfying (cid:104) β n (cid:105) = 0 ; (cid:104) β n β l (cid:105) = 2 αk B T dtδ n,l , (3.10)and the damping coefficients of the algorithm are b = 11 + ( α dt/ m ) , a = b [1 − ( α dt/ m )] . (3.11)The GJF integrator is chosen because of its robustness against discretization time errors, whichis critical for achieving accurate statistics of configurational results. More specifically, it accom-plishes statistical accuracy for configurational sampling of the Boltzmann distribution in closedsystems; and it also provides the correct Einstein diffusion, (cid:104) x (cid:105) = 2( k B T /α ) t , of a freely diffus-ing particle in an unbounded system with constant α [35, 40, 41, 42].We note that Langevin dynamics is diffusive only on time scales larger the so called bal-listic crossover time τ ballistic = m/α , whereas it is predominantly ballistic (inertial) on muchsmaller time scales. Generally speaking, the GJF integrator can be implemented in simula-tions with relatively large time steps, dt > τ ballistic , and still produce accurate statistical resultsat asymptotically large times [35]. A criterion for choosing dt can be set by the requirementthat the characteristic variations in f ( x ) , during the time step, should not be significant, i.e., | f n +1 − f n | (cid:28) | f n + f n +1 | / . This criterion becomes meaningless when a KK interface iscrossed because the interface exerts a singular, delta-function, force [31]. Nevertheless, we willdemonstrate that an accurate algorithm can be devised provided that the integration is performedin the inertial regime with dt (cid:28) τ ballistic (see next section). This implies that the integration timestep in multi-layer systems is bounded by the ballistic time at the most viscous medium: dt (cid:28) τ minballistic = m max( α i ) = min( D i ) mk B T . (3.12)
Before presenting the algorithm for crossing a KK type IBC, the following quantities must beintroduced: • The thermal velocity of the particle, which is independent of α , is given by v th = 2 (cid:90) ∞ vρ MB ( v ) dx = (cid:114) k B Tπm , (3.13)where ρ MB ( v ) is the equilibrium Maxwell-Boltzmann velocity distribution (3.7). • The crossing probability is related to the membrane permeability P and to the thermal ve-locity v th by [31] p = 2 P P + v th (3.14)6 At the interface we have a step-function chemical potential φ step = k B T ln( σ ) H ( x ) , (3.15)where H ( x ) = x < x > (3.16)is the Heaviside step function. The step-function potential result in a delta-function force f step = − dφ step /dx = − k B T ln( σ ) δ ( x ) with a singularity at the interface. In the proposedcomputational scheme, the singular delta-function force is replaced with a sharp, piecewiseconstant force f ( x ) = − k B T ln( σ )2∆ − ∆ < x < − k B T ln( σ )2∆ < x < ∆ , (3.17)defined in the ”small” interval [ − ∆ , ∆ ] = (cid:20) − γD v th , γD v th (cid:21) , (3.18)with the associated potential φ ( x ) = − (cid:90) x −∞ f ( y ) dy. (3.19)The thickness of interface layer (IL) [ − ∆ , ∆ ] over which the chemical potential changesby k B T ln( σ ) is controlled by the dimensionless parameter γ . In the simulations, γ is takento be of the order of unity such that ∆ i ( i = 1 , ) is comparable or smaller of the particlemean free path, l MFP = 2 D i /v th , i.e. the characteristic distance traveled by the particlewithin the ballistic time τ ballistic . The condition (3.12) guarantees that the discrete-timetrajectory does not hop from side to side of the interface, but rather passes across the IL andexperiences the influence of the force (3.17). • We define the weight function W ( x ) = exp (cid:20) φ ( x ) − φ step ( x ) k B T (cid:21) . (3.20)One can easily check that W ( x ) = 1 when f ( x ) = 0 .7 tart a new trajectory (𝑥 ,𝑣 ) Advance trajectory by one time-step 𝑑𝑡 with 𝛼(𝑥 𝑛 ): (𝑥 𝑛 ,𝑣 𝑛 ) → (𝑥 𝑛+1 ,𝑣 𝑛+1 ) No Are 𝑥 𝑛 and 𝑥 𝑛+1 located on different sides of an interface? yes Draw a random number 𝑅 ∼ 𝑈(0,1)
𝑅 < 𝑝 ? [Eq. (3.14)] No Reflect the particle backward (𝑥 𝑛+1 ,𝑣 𝑛+1 ) → (−𝑥 𝑛+1 ,−𝑣 𝑛+1 ) yes Calculate the average friction of the ballistic trajectory, 𝛼 𝑒𝑓𝑓 [Eq. (3.21)] Advance trajectory by one time-step 𝑑𝑡 with 𝛼 𝑒𝑓𝑓 : (𝑥 𝑛 ,𝑣 𝑛 ) → (𝑥 𝑛+1 ,𝑣 𝑛+1 ) 𝑡 𝑛+1 = 𝑡? No Yes Weight by
𝑊(𝑥 𝑛+1 ) [Eq.(3.20)] Update the histogram
Yes Generate another trajectory? No Normalize the histogram Stop Figure 2: A flowchart representation of the algorithm for Langevin dynamics simulations of a layered system with aKK interface.
With the above in mind, the algorithm for calculating c ( x, t ) proceeds as follows:1. Start a new trajectory. Set t = 0 and n = 0 . Choose the initial coordinate x from the initialdistribution c ( x, , and the initial velocity v from the equilibrium Maxwell-Boltzmannvelocity distribution (3.7).2. Advance the trajectory from ( x n , v n ) to ( x n +1 , v n +1 ) by one step dt according to Eqs. (3.8)-(3.11), with f ( x ) given by Eq. (3.17). Use the friction coefficient α ( x n ) at x = x n .3. If x n and x n +1 are found on different sides of the interface then x n +1 needs to be recomputedas follows: We exclude the limit cases σ = 0 and σ → ∞ , which correspond to a perfectly absorbing IBC. The transitionacross such an interface is handled differently, see section 4. Choose a random number, R , uniformly distributed between 0 and 1. • If R > p [with p given by Eq. (3.14)], reflect the particle back to the layer from whichit arrived and set ( x n +1 , v n +1 ) → ( − x n +1 , − v n +1 ) • If R < p , allow the particle to move to the adjacent layer, and determines x n +1 asfollows:3.1 Calculate the ballistic position x n +1 b = x n + v n dt α eff = α ( x n ) | x n | + α (cid:0) x n +1 b (cid:1) (cid:12)(cid:12) x n +1 b (cid:12)(cid:12) | x n | + (cid:12)(cid:12) x n +1 b (cid:12)(cid:12) (3.21)3.3 Advance the trajectory from ( x n , v n ) to ( x n +1 , v n +1 ) by one step dt according toEqs. (3.8)-(3.11), with the effective friction coefficient α eff (3.21). Notice thatin some rare cases, the new position x n +1 will be found on the same side as x n ,but this is acceptable since small discretization errors are always present whenencountering a step function diffusion function.4. If t n +1 = t then • Stop the trajectory at x = x n +1 . • Weight it with the weight function W ( x ) (3.20), and update the histogram , hist w ( x ) ,for the distribution function w ( x, t ) : hist w ( x ) = hist w ( x ) + W ( x ) . • Return to step 1 if you want to generate another trajectory; otherwise go to step 6.5. Return to step 2.6. Normalize the distribution, w ( x, t ) , to obtain the concentration profile, c ( x, t ) : c ( x, t ) = w ( x, t ) (cid:82) ∞−∞ w ( x, t ) dx . (3.22)Figure 2 shows a summary of the algorithm in the form of a flowchart. In this section we consider a biomedical example where the previous concepts and algorithms areapplied to a simple model of a drug-eluting stent (DES). Stents are small mesh tubes inserted tokeep open stenosed arteries (see fig. 3). Drug-eluting stents (DES) also have an additional thinlayer of polymer coating the mesh and eluting a drug. More precisely, a DES is constituted by In the histogram representation, hist w ( x ) , data accumulate in discrete bins. The continuous distribution w ( x, t ) is defined as the total value stored within the relevant bin, divided by the bin size. igure 3: A drug-eluting stent implanted in an artery. metallic prosthesis ( strut ) implanted into the arterial wall and coated with a thin layer of bio-compatible polymer that encapsulates a therapeutic drug ( coating ). Such a drug, released in acontrolled manner through a permeable membrane ( topcoat ), is aimed at healing the vascular tis-sues or at preventing a possible restenosis by virtue of its anti-proliferative action against smoothmuscle cells [43, 44].To formulate the mathematical problem that serves as a simple DES model, let us consider astent coated by a thin layer (of thickness L ) of polymer containing a drug and embedded intothe arterial wall (of thickness L ), as illustrated in fig. 4. The complex multi-layered structure ofthe arterial wall has been disregarded for simplicity, and a homogeneous material with averageddiffusion coefficient D has been considered. A small plasma filtration velocity is present in thewall, but a scaling analysis shows that this transport effect remains negligible in comparison withthe diffusive one [43, 47]. The diffusion coefficient of the polymer is D (cid:28) D . The DES modelshown schematically in fig. 4 is a two-layer system similar to the one depicted in fig. 1. Theonly difference between them is that here the two layers have a finite extent and two boundaryconditions (BCs) are prescribed to make the mathematical problem well-posed. Since the strutis impermeable, no mass flux passes through the left boundary surface, which is modelled byimposing a reflecting boundary condition: J ( − L ) = 0 . The right side L , being L (cid:29) L , ismodeled as an absorbing boundary, namely c ( L ) = 0 . At the initial time ( t = 0 ), the drug iscontained only in the coating (layer 1) and it is uniformly distributed at a maximum concentration C : c ( x,
0) = C for − L ≤ x ≤ c ( x,
0) = 0 for ≤ x ≤ L . (4.1)To slow down the drug release rate, a thin membrane (called topcoat ) is located at the interface x = 0 between the two layers. The topcoat separating the coating and the arterial wall imposesthe KK IBC (3.6) between the layers. As no drug is lost in the topcoat, the continuity of the fluxIBC (3.5) is also assumed.To summarize, the two-layer diffusion problem is given by the following set of partial differ-10 (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (a) (b) (c) (d) c x L c (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) lumenwall (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) Figure 4: Cross-section of a stented artery with the sequence of layers for drug dynamics (a) stent strut, (b) coating,(c) topcoat, (d) arterial wall (figure not to scale). ential equations, with boundary and initial conditions [47]: ∂c ∂t − D ∂ c ∂x = 0 in [ − L , (4.2) ∂c ∂t − D ∂ c ∂x = 0 in [0 , L ] (4.3) − D ∂c ∂x = − D ∂c ∂x = P ( c − σc ) at x = 0 (4.4) ∂c ∂x = 0 at x = − L (4.5) c = 0 at x = L (4.6) c = C, c = 0 at t = 0 (4.7)The solution of the above problem is obtained by separation of variables: c i ( x, t ) = X i ( x ) G i ( t ) i = 1 , (4.8)11here the spatial functions X and X satisfy the Sturm-Liouville problem: X (cid:48)(cid:48) = − λ X in [ − L , (4.9) X (cid:48) = 0 at x = − L (4.10) D X (cid:48) = D X (cid:48) at x = 0 (4.11) X (cid:48)(cid:48) = − λ X in [0 , L ] (4.12) X = 0 at x = L (4.13) − D X (cid:48) + P σX = P X at x = 0 (4.14)with: D λ = D λ (4.15)The general solution of the ordinary differential eqns. (4.9) and (4.12) is: X ( x ) = a cos( λ x ) + b sin( λ x ) X ( x ) = a cos( λ x ) + b sin( λ x ) (4.16)and G ( t ) = G ( t ) = exp( − D λ t ) = exp( − D λ t ) (4.17)The eigenvalues λ i and the unknown coefficients a i and b i are computed by imposing the BCs andIBCs as follows. From (4.10) and (4.13), we have: a sin( λ L ) + b cos( λ L ) = 0 a cos( λ L ) + b sin( λ L ) = 0 , (4.18)and from (4.11) and (4.14), it follows that D b λ = D b λ − D b λ + P σa = P a . (4.19)Eq. (4.18)–(4.19) form a system of four homogeneous linear algebraic equations in the four un-knowns a , b , a and b . To get a non trivial solution, it is needed that the determinant of thecoefficient matrix associated with the above system be equal to zero, that is: tan (cid:32)(cid:114) D D L λ (cid:33) ( D λ + P σ tan( λ L )) − (cid:114) D D P = 0 (4.20)An infinite sequence of eigenvalues λ , λ , ....λ m ... is obtained as solutions of the above tran-scendental equation (4.20) (eigencondition). Hence, the complete solution of the problem (4.2)–(4.7) is expressed as a linear superposition of the fundamental solutions: c ( x, t ) = ∞ (cid:88) m =1 A m X m ( x ) exp( − D λ m t ) c ( x, t ) = ∞ (cid:88) m =1 A m X m ( x ) exp( − D λ m t ) (4.21)12 igure 5: Drug concentration profiles in the coating (above) and in the wall (below) for three times (note the differentspace scales). The curves depict the solution obtained by separation of variables in [47], while the markers representthe results of the Langevin simulations based on the algorithm described in section 3.2. where A m are determined through the initial conditions (4.7) (see [47] for further details). In the absence of direct experiments, we have chosen the following parameters which are in thecorrect range and for which the resulting release times are consistent with published data [44, 45,46]: L = 5 · − cm L = 10 − cm D = 10 − cm /s D = 7 · − cm /sP = 10 − cm/s σ = 0 . (5.22)These parameters, which are representative of the typical scales in DES, have been chosen basedon data in literature for the arterial wall and heparin drug in the coating layer. The same parameterswere used in ref. [47], with the exception of D and D have been taken smaller, in order to13 igure 6: The concentration c at t = 10 in the region close to the interface. The markers depict the Langevinsimulation results with different value of γ in Eq. (3.18), while the dashed red line presents the analytical solution inthe limit γ → , when the chemical potential is given by a step-function - see Eq. (3.15). have more realistic release times. For the Langevin simulations, we use dimensionless units with k B T = 1 , v th = 2 . , m = 10 − , L = 5 , L = 100 , D = 10 − , D = 7 . For γ = 0 . inEq. (3.18), ∆ (cid:39) · − (cid:28) L and ∆ = 1 . (cid:28) L . In these units, P = 10 − . Converting thedimensionless units to physical ones, we find that t = 1 in the simulations corresponds to s .The time step is set to · − , which falls in the ballistic regime of the Langevin dynamics inboth layers, τ minballistic = 10 − [see Eq. (3.12)]. We note that the reflecting boundary at x = − L istreated as special cases of the KK condition with P = 0 and σ = 1 and is, therefore, covered bythe above algorithm. The absorbing boundary at x = L corresponds to P → ∞ and σ → ∞ (or σ = 0 ). In this case, one should assign a very large (or nearly vanishing) value for σ in Eq. (3.17).In our simulations we use a simpler approach: We do not introduce a force near the absorbinginterface and, instead, simply terminate and assign zero weight to each trajectory exceeding L .The concentration profiles, c and c , for three values of time are displayed in fig. 5. Weobserve that the concentration c decays in time, indicating that drug is eluting from coatingto the wall. The concentration at the wall, c , increases at short times, and decays at longertimes as more and more drug arrives at the absorbing surface x = L . At x = 0 , wherethe KK IBC is imposed, we observe a sharp discontinuity in the concentration that diminisheswith time. The agreement between the semi-analytical solution (continuous curves) and theLangevin simulation results (diamond symbols) is excellent, except for deviations near x = 0 at the shorter time t = 10 s . These arise from the approximation of the delta-function forceat the KK interface by the sharp continuous force (3.17) existing around the interface. Theimpact of this approximation on the results are supposedly corrected by the weight function143.20); however, this correction is based on the ratio of the corresponding Boltzmann factorsand, thus, relies on the assumption that locally the system is at thermal equilibrium which, strictlyspeaking, can be only assumed in the overdamped limit τ ballistic → . Fig. 6 presents resultsfor c at t = 10 with larger values of γ , zooming in on the region close to the interface.The difference between the analytical and numerical solution at x = 0 provides a measure ofthe computational error, Er . Not surprisingly, we find that it decreases almost linearly with γ ( γ = 4 : Er = 0 . , γ = 2 : Er = 0 . , γ = 1 : Er = 0 . , γ = 0 . Er = 0 . )suggesting that the simulations should be run with the smallest possible γ . Nevertheless, γ cannotbe reduced indefinitely since the condition dt (cid:28) γl MFP / v th is required to ensure that the particletravels within the IL . In the last section, we have examined the computational error arising from the approximationof a discontinuous chemical potential with a sharp piecewise constant jump. Here, we furtherexpand our analysis, focusing on the convergence and accuracy of the algorithm with respect tothe integration time step dt . As noted above, we use the GJF equations (3.8)-(3.11) to integratethe Langevin dynamics, where an ensemble of particle starting on one side of the interface andspreading across the system. We chose this integrator because it yields the correct Einstein dif-fusion, (cid:104) x (cid:105) = 2 D t = 2( k B T /α ) t , for any time step when applied in simulations of a freelydiffusing particle. Thus, the algorithm samples correctly the diffusive dynamics away from theinterface, and discretization errors arise from the segments of the trajectories when the particlepasses close to the interface. These errors can be minimized by using smaller dt , but that wouldcome at the cost of being able to simulate a smaller number of trajectories per CPU time, whichwould increase the statistical noise. In order to analyze the convergence of the numerical methodwith respect to dt , we repeat the simulations of a system with IL parameter γ = 0 . for a sequenceof decreasing time steps dt k ( k = 0 , , , .... ). As a reference case, we set dt = 25 · − whichis 50 times larger than the minimal time step dt used to generate the results in fig. 5 and 2.5 timeslarger than the ballistic time, as computed from Eq. (3.12). We quantify the distance between theconcentration profiles c ( k ) corresponding to subsequent time-steps through the Euclidean norm E k ( · , t ) = (cid:13)(cid:13)(cid:13)(cid:13) c ( k ) − c ( k − c ( k ) (cid:13)(cid:13)(cid:13)(cid:13) , k = 1 , . . . , (6.23)The results of the analysis are summarized in table 1. The table shows a clear convergence atsmaller time steps and indicates that choosing dt = 5 · − for the simulation results in fig. 5yields a satisfactory accurate solution. The significant drop in E k between k = 2 and k = 3 isprobably due to the fact that dt is not sufficiently smaller compared to the ballistic time ( dt =5 · − = τ ballistic / ). Thus, for the smaller k values in the table the error is predominantly asystematic discretization one, while for the larger values of k is dominated by statistical noise. Note that the above condition can be also written as dt (cid:28) γ ( π/ τ ballistic , with τ ballistic given by Eq. (3.12),which explains why γ should be of the order of unity. E ( · , ) E ( · , ) E ( · , )1 dt = 20 dt dt = 10 dt dt = 5 dt dt = 2 dt dt = dt Table 1: Norm of the profile difference at three times for a sequence of decreasing time steps dt k . To summarize, the simulation results shown in fig. 5 represents an acceptable compromisebetween accuracy and computational efficiency, dictated by the available CPU time, the highaspect ratio ( L /L = 20 ), and the large diffusivity contrast ( D /D = 700 ). We proposed an algorithm for Langevin dynamics simulations in diffusive multi-layer systems,with flux continuity and KK interface condition separating regions of different diffusivity. Theproposed method is based on accumulating statistics from a large number of independent singleparticle trajectories. These are produced by a Langevin dynamics discrete-time integrator, and theproposed algorithm describes how the integration is set up when the particle crosses an interface.From the ensemble of Langevin dynamics trajectories, we generate a fine-grained histogram ofthe concentration profile that solves the corresponding continuum diffusion equation.To validate the algorithm, we consider the case study of two-layer model for a DES that can besolved semi-analytically by separation of variables. The agreement between this solution and ourcomputational results is shown to be very good. We also use this example to assess the accuracyand stability of the method. Our analysis suggests that two parameters of the simulations need tobe carefully chosen: (i) The integration time step that must be smaller than the ballistic time of theLangevin dynamics, and (ii) the width of the interface layer over which the step-function potentialenergy is approximated. Reducing the values of these parameters improves the accuracy of theresults, but also increases the computational cost since more iterations are needed for generatingeach trajectory. A careful choice should balance between these two aspects, and depends on theproblem in question and on the available computational resources.While the example discussed here concerns a two-layer system, it should be stressed that aclear advantage of the Langevin dynamics algorithm is in dealing with multi-layer systems thathave relevance to applications in many scientific and engineering disciplines. The method canbe straightforwardly generalized to any number of interfaces, simply by employing the algorithmwhenever a trajectory encounters one of the interfaces. The simplicity of the algorithm is incontrast to analytical solutions that, in general, become increasingly complex and computationallyinefficient with larger number of layers. In a future work we plan to present studies of multi-layered systems to demonstrate this important feature of the method.Another direction is to extend the method to two- and three-dimensional composite systems.16e also intend to consider examples where other mechanisms besides passive diffusion, e.g. ad-vection and mass degradation, are included. For the specific application of drug-eluting stentconsidered herein, additional efforts are needed to assess and evaluating the relative influence ofthe various factors, including material properties.
Acknowledgments
Funding from the European Research Council under the European Unions Horizon 2020 Frame-work Programme (No. FP/2014-2020)/ERC Grant Agreement No. 739964 (COPMAT) is ac-knowledged.