Nanomedicine by extended non-equilibrium thermodynamics: cell membrane diffusion and scaffold medication release

1.   Introduction

Nanomedicine knows an increasing interest due to its estimated importance in efficient medication delivery ^[1]. Crucial medical topics deal with its application in several fields, such as overcoming the blood-brain barrier ^[2], anti-diabetic medication ^[3] and uptake in endothelial cells ^[4]. Nanomedicine deals with the novel chemical, physical or biological properties of medication at nanometric scale, as well as the application of nanotechnology aiming at making a medical diagnosis, treating or preventing the cause of illness. Therefore, nanomedicines have the ability, due to their particular size or surface functionality, to demonstrate physicochemical and/or biological properties that make them interesting for medication delivery, controlled medication release and/or improved medication transport through biological barriers, such as cell membranes, which would otherwise be much more difficult ^[5]. One of the key parameters for successful medication delivery is controlling or improving the permeation of medication through the destined place (e.g. vascular system, cell) in the human body ^[6,7,8]. When speaking about cell permeation and the body uptake efficiency of drugs into the body that it entails, many works take the apparent cell permeability, P_app, through monolayers of human epithelial cells as a widely accepted standard ^[9,10,11]. Values of P_app (often presented in medicine as a volumetric flux per unity of area) of such cells are also often used for the estimation of the blood-brain barrier ^[12,13]. It is useful to understand the apparent cell permeation as a series of parallel and serial resistances (from layer to another). This vision is based on the concept of solubility diffusion ^[14]. The permeability is seen here, therefore, not as an imposed force, but is rather passive in a way that only spontaneous diffusion governs the permeation process. This diffusion process is determined by a concentration gradient. The concentration gradient is classically described by Fick's law or, in case of sizes (within the membranes separating the monolayers) at nanoscale, by a more general evolution equation, describing the extended diffusion law, as will be discussed later. Active transport, as could be caused by transporter proteins, is not considered here. Although, in reality, chemicals/drugs could be transported by proteins (and thus have an effect on P_app), it is not taken into account here. It is well admitted that generally three pathways are possible through a cellular monolayer, i.e. in between the cells (the paracellular way), through the cell (considered here) and along the membrane solely (the lateral way) ^[15,16]. In this case, parallel resistances of the three pathways would apply. Here, we only consider the pathway through the cell. Following the pH partition hypothesis ^[17], it can be assumed that only the neutral form of the drug chemicals can pass the membrane resistances. Therefore, chemicals that are ionized should undergo an acid-base reaction in between the water layer and the membrane before they are able to pass the membrane as a neutral species. According to the aforementioned pH partition hypothesis, such a reaction can be considered to occur instantaneously, meaning that no kinetic hinderness (and therefore no additional resistance) needs to be considered.

It is often suggested that the diffusion through membranes are of non-Fickian nature. We propose that this is either because of non-local effects (the size of the membranes are of the same order of magnitude as the drug molecules or not much bigger) or because the relaxation time becomes important. As controlling and enhancing drug delivery is an important issue in the field of nanomedicine, next to in vitro experiments, modelling is crucial to understand the underlying mechanisms. In nanomedicine, the inertia of the diffusing particles is relevant in addition to the small size of the channels through which the medication passes, being in some cases of the same order of magnitude as the mean free path of the particles. Furthermore, purely microscopic models are often difficult to implement directly to real-life applications. Therefore, further progress in this field can be enhanced by models that are easily applicable, whilst keeping a qualitative efficient and understandable way to link to experiments in a multi-disciplinary way. Therefore, we propose here a model based on thermodynamic principles, where we extend them by incorporating effects, in a rigourous way, that account for the relaxation and size of the drug molecules with respect to the channels through which they diffuse. With the purpose to demonstrate the applicability of our new theory, we treat two cases: drug diffusion through a cell membrane and drug release out of a scaffold. For the first case, let us imagine the existence of two compartments of well-mixed water that are separated by a barrier, i.e. a membrane: the first containing the drug molecules as the donor and the second represents the inner cell or the receptor. The purpose is to follow the evolution of the mass fraction profile through the membrane at the beginning of the diffusion process, where we assume that the drug molecules (in that beginning process) concentration is that high that the mass fraction at the inlet of the membrane is constant. For the second case, we discuss drug release out of a scaffold and follow the concentration profile to assess the rapidity of the drug release. To account for the importance of the relaxation time of the drug molecules and the non-local effects, we propose the formulation of mass transport using Extended Non-Equilibrium Thermodynamics ^[18,19], whose main characteristic is to upgrade the flux and higher order fluxes to the rank of independent variables.

2.   Extended equation for the mass flux through a cell membrane

At micro and nanoscales, mass transport is mostly influenced by non-local effects and relaxation times. The classical Fick law

relating the mass flux vector J to the mass fraction gradient ∇c, with ρ and D denoting the density and mass diffusion coefficient, respectively, is not applicable at short diffusional time (with respect to the relaxation time, i.e. large relaxation times) and small spatial scales. In order to account for large relaxation times, Fick's law can be generalized by a Cattaneo-like equation (in analogy to the one used for heat transport ^[20]) under the form

where τ_r designates the relaxation time of the mass flux and ∂_t the partial time derivative. To account furthermore for non-local effects, Extended Non-Equilibrium Thermodynamics (see for underlying developments in previous work ^[18]) proposes a further generalization

It is convenient to reformulate (3) in dimensionless form. Therefore, the mass flux J is rescaled with $\frac{{\rho D}}{\delta }$, the space coordinate with a characteristic length δ (standing also for the membrane thickness), the time t with the characteristic diffusion time ${{\rm{ \mathsf{ τ} }}_d} = \frac{{{\delta ^2}}}{D}$. Introducing the Knudsen number $Kn \equiv \frac{\ell }{\delta }$, a dimensional number $M = \frac{{{\tau _r}}}{{{\tau _d}}}$ and using the same notation for the non-dimensional quantities, Eq (3) becomes

It should be noted here that by choosing τ_d as time scale, one implies that any change in M means a change in the relaxation time, which evokes non-local effects. This motivates proposing an adaptation for the expression for M. Indeed, kinetically speaking, one may define a relation $D \equiv \frac{1}{3}\frac{{{\ell ^2}}}{{{\tau _r}}}$, in analogy to the thermal conductivity ^[1]. Using the definition of τ_d and Kn, one finds a relation $M = \frac{1}{3}K{n^2}$. However, in suspensions, these kinetic considerations are not always valid and differences could occur. Acknowledging that some partial proportionality is admitted between M and the Kn number (the relaxation time and the mean free path are related in that they both can cause non-local effects), this motivates to define rather $M = {f_{Kn}}\frac{1}{3}K{n^2}$. Note that using this definition does not change anything for the analysis to be performed (indeed, in principle, any number could be chosen for M in the framework of a dimensionless analysis conditioned by correct interpretation of the results). Moreover, f_Kn is a factor that accounts for the observation that in nanofluids the relaxation time is not necessarily fully proportional to the Kn number. This factor also allows investigating separately the effects of the relaxation time and the mean free path on drug diffusion.

We can see clearly differences between (3) and the conventional constitutional Fick's law for mass diffusion (1), but also some similarities. When we deal with systems where the characteristic size is found to be much larger than the mean free path of the molecules l, Eq (3) reduces to the Cattaneo-like equation (2). This would also correspond to neglect any higher order of mass fluxes. Now, when the changes in the mass flux occur in a time span that is much larger than the relaxation time, we obtain Fick's law. Let us stress that the latter issue depends on the physicochemical relation the medication molecule can have with the carrier fluid or the host environment. Other well-known methods that are used for solving models that deal with nanoscale systems are the phonon-Boltzmann equation, Monte Carlo simulations and Molecular Dynamics. The Monte Carlo simulations and Molecular Dynamics method are numerical methods that generally do not require higher order variables, but instead demand high computing times, whilst being less flexible for modelling various scales at a time, although some efforts have been performed concerning multi-scale molecular dynamics ^[21]. Nonetheless, the flexibility in scale remains less than the present approach. The phonon-Boltzmann equation also demands high computational time for solving a model. The aforementioned aspects are generally less of a problem when it comes to the model that is proposed in this work. It should be noted that the present development could be extended by introducing non-linear mass fluxes or higher-order (higher than the second order) fluxes. However, implementing such an extension, albeit theoretically not an impossible task (see ^[18]), would introduce numerical difficulties and higher computing time, which is an absent aspect of the present approach. Furthermore, it would also loose the ability to propose analytical solutions for some particular cases, offering a pedagogical presentation allowing easier interpretations. Non-linear contributions of the mass flux (at first and higher orders) could be interesting to consider when the length scale becomes much smaller than the mean free path of the molecules. For the purposes of nanomedicine applications, this would in many cases lead to unnecessary complications. Indeed, for many cases a linear mass flux limited to second-order mass fluxes (which leads to the term ∇²J in (3)) is sufficient. In comparison, it already appeared to be sufficient for heat transport applications ^[22] and for predicting nanoflow properties ^[23].

3.   Nanomedicine: non-local diffusion

3.1.   Membrane permeation

From typical clinical observations that lead to data in the literature, it appears that one of the most important issues to be studied for improving drug delivery are through enhancing the permeability. In porous media, the concentration gradient can be negligible with respect to the flux, which is often the case for nanoporous flow that is pressure-driven by some external pressure. There are some cases in the human body that are comparable from a conceptual point of view, such as oxygen and carbon dioxide transfer from the lungs to the veins within the alveoles. Nonetheless, porous flow within cell membrane pores presents often the difference that, due to the absence of externally imposed pressure, the concentration gradient can be quite important with respect to the flux. In case of passive permeability, one speaks often of diffusional processes. As the diffusion process through the membrane is often the limiting pathway, the flux across the membrane will undergo non-local effects. Flux variation in the perpendicular direction to the flux across the membrane can be neglected, assuming a well-dispersed medication solution at the entrance of the whole membrane. One would then use (3) to describe such a process, reduced to a one-dimensional equation, replacing the ∇ operator by the derivative ∂_y. However, the concentration gradient is not simply a gradient of the concentrations from both sides of the membrane, since concentration jumps can exist at the entry or exit of those membranes due to lipid bilayer interaction or wettability effects. The gradient that counts is the one within the membrane. Therefore, it is important to consider the fact that flux by diffusion across a membrane (even if it were uniform) depends on the solubility. As there are different solubilities at both sides of the membrane and within the membrane as well, we can define partition coefficients ϑ_in and ϑ_out for the entry and exit of the membrane, respectively:

where N_{mb, in} and N_{mb, out} are the concentrations within the membrane from the entry and exit sides, respectively. Also, N_in and N_out stand for the same concentrations, but at the medium sides (outside the membrane), respectively. It is more convenient to write (5) and (6) in mass fractions, c, where the subscripts mean the same thing as for the concentrations. Assuming ϑ_in = ϑ_out≡ϑ, we have

where ϑ is to be determined experimentally as a material property. We can readily see that for a higher partition coefficient, the concentration gradient within the membrane is then larger than on the outer sides of the membrane, which results into a higher membrane solubility and thus a higher permeation. To model the whole process, next to the one-dimensional version of Eq (3), the mass fraction evolution is given by the classical species balance equation, which in dimensionless form gives

Since the thickness of a cell membrane is much smaller than its size in the other directions, we consider a one-dimensional model through the membrane with thickness δ, assuming that in all directions other than across the membrane the variables remain constant. The entrance of the membrane (at the drug donor side) is denoted by the dimensionless coordinate y = 0 and the exit (at the drug receptor side) by y = 1. The transport of the drug molecules through a cell membrane will then be modeled by the one-dimensional versions of Eqs. (4) and (9)

Here, we omitted the subscript mb for the membrane and will do so in the following. The drug at the donor side of, but outside, the membrane is given by c_in. The inner-side drug mass fraction of the membrane at y = 0 is then, via the partition theory, given by

At the cell-side of the membrane, y = 1, we assume a continuity of fluxes so that

where the subscript ce denotes the inner cell. Assuming that the inner cell is an unstirred well-mixed environment, the drug uptake can be considered to occur under a constant rising concentration. Therefore, the gradient is nearly zero, so that we can say that $\frac{{\partial {c_{ce}}}}{{\partial y}} \approx 0$. This allows to state approximately that at y = 1

The initial mass fraction within the membrane is given by

3.2.   Medication delivery through a scaffold

A supplementary interest is the continuous delivery of drugs out of a scaffold. In this case, a scaffold, containing the drugs, is introduced into the area of interest (the manner in which this is done is not the point of discussion here). This scaffold is a porous structure, where we are interested to provide a tool for analysing the diffusional flow through this medium. One side of the scaffold contains the drugs, maintained at a certain concentration or mass fraction, c_in. The other side is where the drugs enter the area of interest, where we assume that the uptake is immediate so that the mass fraction is zero there, i.e. c₀ = 0. Let us approximate such a scaffold as cylindrical pores. Contrary to the previous subsection, the dimension of the flux direction (axial direction of the pores) is larger than the radius of such a cylinder (radial direction of the pores). In that case, we can say (and also motivated by previous work in porous systems ^[23]) that the flux in the axial direction of the pores (the direction of drug release) will be approximately constant. However, the flux will be variable over the radial direction of the flux. Indeed, due to the small pore size with respect to the axial dimension of the pores, the wall effects will influence the flux and a flux profile exists in the radial direction. As for the mass fraction gradient (being a slow variable), its radial gradient is negligible with respect to the axial one (influenced by the drug release). As such, Eq (4) gives in dimensionless form (using $\frac{{\rho D}}{R}$, R and L as scaling for the mass flux, the radial spatial coordinate and the axial spatial coordinate, respectively, where R is the pore radius and L the length of the pore)

with $K{n_r} \equiv \frac{\ell }{R}$, r the radial spatial coordinate, x the axial spatial coordinate in the pores and ${\rm{ \mathsf{ α} }} \equiv \frac{R}{L}$. Assuming maximum flux in the middle of the pore (r = 0), i.e. ${\left. {\frac{{\partial \mathit{\boldsymbol{J}}}}{{\partial r}}} \right|_{r = 0}} = 0$, and first-order slip (defining ${B_s} \equiv \frac{{{\ell _s}}}{R}$, where ℓ_s stands for the slip length, related to the wall material) at the wall (r = 1), i.e.

the flux can be found as a function of r. Taking subsequently an average value over a cross-section of the cylinder gives the following mass flux

with ${\cal B}$ the Bessel-I function. It should be noted that Eq (18) depends on how slip is observed and how the mass flux is defined. Other effects can influence the mass flux, but here it is the purpose to show how physical models can be used in the field of nanomedicine. Such an effective flux (taking into account the material of the scaffold (via the slip factor) and non-local effects via the Knudsen number) is now filled in the dimensionless (time scale is taken as $\frac{{{L^2}}}{D}$) species balance in the flow direction

The evolution of the mass fraction is then described as

The boundary condition on one side of the scaffold is given by

at x = 0, which stands for assuming that the drug only goes on one direction and

at x = 1, which allows for a drug uptake relaxation time of the same order of magnitude as the characteristic diffusion time (which is reasonable, considering that transport in cells is mainly of diffusive character), describing the mass fraction going from c₀ to 0 (complete drug uptake). This depends of course of the system in question and the area in the body. Nevertheless, studying the effect of non-locality and slip factor on drug uptake, the choice of a relaxation time is not of importance here and does not influence our analysis. Eq (22) also shows the initial condition by evaluating it at t = 0.

3.3.   Numerical method

The model is solved using a semi-implicit difference method for the spatial discretization and the forward Euler method for the temporal discretization for each time step t. At the beginning of the simulation, the initial conditions are imposed. At a certain time t, the equations for the mass flux and mass fractions are written in matrix form A*V = B, where A is a matrix of dimension (n-2)×(n-2) (n being the number of nodes in the discretized spatial domain and two nodes being dedicated to the boundary conditions, the bulk equations are thus discretized on (n-2) nodes) containing the coefficients in the equations for the mass fractions and mass fluxes at time t+1, V a vector of dimension (n-2) containing the unknown mass fraction and mass flux values at the spatial nodes at time t+1 and B a vector of dimension (n-2) containing the known mass fraction and mass flulx values at time t. Due to the spatial gradients in the bulk equations, the values at nodes 2 and n-1 depend on the boundary values at nodes 1 and n, respectively. The vector V is calculated as V = A^-1B. The boundary condition for at the membrane entrance is imposed and held at that value at each time step. The boundary condition at the exit of the membrane at time t+1 is then obtained through the bulk values at time t+1. This procedure is repeated until the relative errors of the bulk and boundary values of both the mass fraction and mass flux between the previous and present loop become smaller than 10^-4. The obtained values at time t+1 are stocked in the output matrix and used as the known values for the next time step t+2.

4.   Results and discussion

4.1.   Membrane diffusion

We use our model, given by Eqs. (10)–(12), (14)–(15), and solve it by means of a semi-implicite finite difference algorithm, explained in section 4.1. We take c_in = 0.8, c₀ = 0 and ϑ = 1 (we neglect the effect of the partition coefficient for the moment). The results are presented in Figures 1-5. Figures 1-3 show the evolution of the mass fraction across the membrane for Kn numbers Kn = 0.1, Kn = 1 and Kn = 10, respectively. The different Kn numbers stand for different sizes of the membrane with respect to the molecular mean free path of the drug molecules in the fluid present in the membrane (which can be taken as a watery solution). Furthermore, for each Kn number different f_Kn values are used: f_Kn = 0.1, f_Kn = 1 and f_Kn = 10.

Figure 1.  Mass fraction evolution across the membrane layer for different values of f_Kn for Kn = 0.1.

DownLoad: Full-Size Img PowerPoint

Figure 2.  Mass fraction evolution across the membrane layer for different values of f_Kn for Kn = 1.

DownLoad: Full-Size Img PowerPoint

Figure 3.  Mass fraction evolution across the membrane layer for different values of f_Kn for Kn = 10.

DownLoad: Full-Size Img PowerPoint

Figure 4.  Mass fraction as a function of time at y = 1 for different values of f_Kn and Kn = 0.1 (left), 1 (middle) and 10 (right).

DownLoad: Full-Size Img PowerPoint

Figure 5.  A zoom around y = 1 for several f_Kn values for Kn = 1 (left) and for several Kn values for f_Kn = 1 (right) at y = 0.3.

DownLoad: Full-Size Img PowerPoint

We can see that for a given Kn number a higher f_Kn number (higher relaxation time with respect to a fixed diffusion characteristic time), the evolution of the mass flux is significantly slower, with a large gradient at the membrane entrance. This implies that the mass flux is very high at the membrane entrance but very low at the membrane exit. This effect is precisely referred to as non-locality. It is important here that by choosing the diffusion characteristic time as the scale time, any change in f_Kn implies a change in the relaxation time. This means that one compares systems with different relaxation times, which illustrates the effect of medication diffusion as a function of different medication-carrier fluid systems. A higher Kn number has the same non-local effect, albeit on the geometrical level (for the same system). This latter effect is not so clearly visible in Figures 1–3, but can be more appreciated looking to the evolution of the mass fraction at the membrane exit y = 1 in Figure 4.

Effectively, Figure 4 shows that a higher Kn number, for all f_Kn values, results into a lower mass fraction, which is representative for a lower drug diffusion/solubility in the membrane. The non-local effect at high Kn numbers is responsible for a non-Fickian mass fraction profile, which is indicated by a mass fraction at the exit of the membrane that is less than the initial mass fraction at steady state (this in comparison with Fickian behaviour, where for a zero-flux boundary condition, the profile should become flat). This means that for drugs that have large relaxation times or mean free paths, the drug uptake could be troublesome. It is interesting to note that the values shown in Figure 4 are exactly at the membrane exit, but that a mass fraction jump exists just before that exit. Not being well discernible in Figures 1-3, a zoom of some of those curves has been plotted in Figure 5.

This mass fraction jump can be seen clearly in Figure 5 and stems from the non-local behaviour from both the mass flux time derivative (∂_tJ) and mass flux diffusion term (∇²J). Now, if the partition coefficient was not at unity, we would have not only a jump due to a non-locality at the membrane side of the interface, but also at the cell side, due the partition coefficient. Figure 6 shows the influence of the partition coefficient on the mass fraction evolution.

Figure 6.  Mass fraction across the membrane at t = 3 for f_Kn = Kn = 1 (left) and as function of time at y = 1 for f_Kn = Kn = 1 (right).

DownLoad: Full-Size Img PowerPoint

It can be seen that a higher partition coefficient results into a higher mass fraction and a higher drug uptake (higher mass fraction at the membrane exit). Since a zero-flux boundary condition (standing for a quick drug diffusion in the cell and not for a wall-like boundary, as is discussed under Eq (13)) is maintained at the membrane exit, the partition coefficient does not affect the mass fraction here. However, it does effect the quantitative value of the mass fraction just outside the membrane (in the cell). The higher the partition coefficient, the lower the mass fraction will be here. Also, as Figure 6 seems to show, a higher partition coefficient does not exhibit a different behaviour of the mass fraction at the membrane exit. That is, a higher partition coefficient gives a higher mass fraction of the drug at the membrane entrance, which translates itself into a higher mass fraction at the membrane exit. Nonetheless, this increase is cancelled by the partition coefficient itself at the membrane exit. This motivates using drugs that have a higher partition coefficient at the membrane entrance, but a lower one at the membrane exit. This would, of course, not lead to a higher total uptake (mass conservation implies that the total quantity of the uptake is equal to the administered quantity), but may lead to a higher efficiency if a higher mass fraction is needed for a better healing in a shorter time period. Having discussed drug uptake by means of drug diffusion through a cell membrane, it is also interesting to provide tools to drug release by scaffolds, applicable in the field of nanomedicine. This is treated in the next subsection.

4.2.   Scaffold drug release

We use our model, given by Eqs (20)–(22), and solve it in the same way as for the non-local drug diffusion model. We take c₀ = 0.8 and present the results in Figures 7 and 8 for different Kn numbers and slip factors B_s. When we take B_s = 0, this stands for a zero slip length of the mass flux. This means that the mass flux at the wall of the pores in the scaffold is equal to zero. However, when B_s→∞, this means that there is full slip at the pore walls and that the mass flux is everywhere equal in the radial direction. We take a value B_s = 10 to be quite large, standing for "large slip" and compare its results to "no-slip". The slip behaviour of the mass flux is not clearly defined in the literature, but is rather often used in the framework of fluid flow, related to its velocity. The very definition of a mass flux in a porous medium being ^[23]

Figure 7.  Mass fraction across the scaffold pores for several times t for Kn = 0.1 (left), 1 (middle) and 10 (right) with a no-slip condition at the pore walls, B_s = 0.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Mass fraction across the scaffold pores for several times t for Kn = 0.1 (left), 1 (middle) and 10 (right) with a high-slip condition at the pore walls, B_s = 10.

DownLoad: Full-Size Img PowerPoint

where ϕ is the porosity of the porous medium, ρ_f the density of the fluid flowing through it, v_f the fluid velocity and v the barycentric velocity, we see that the mass flux is directly related to the fluid velocity. If we note that in a two-component system the barycentric velocity is expressed as ρv = ϕρ_fv_f+(1-ϕ)ρ_sv_s, wherein the total mass density ρ is ρ = ϕρ_f+(1-ϕ)ρ_s with ρ_s and v_s the density and velocity of the solid component, respectively, and that in case of a scaffold v_s = 0, we can quickly see that in (23) the only variable the mass flux is related directly to is the fluid velocity v_f. As such, the discussion on the slip factor becomes more physical and understandable.

Figures 7 and 8 show that the Kn number tends to delay the drug release. This suggests that for pores of sizes of the order of magnitude of the mean free path of the drug molecules the drug release is slowed down. On the other hand, in some applications, it is necessary to have pores of such a size. Figure 8 shows that the slower drug release can be compensated by a higher slip factor, which accelerates the release of medication, illustrated by mass fraction profiles that evolve faster with respect to Figure 7. Figure 9 shows the mass fraction evolution for different Kn numbers for "no-slip" (B_s = 0) and "high-slip" (B_s = 10) at the membrane exit (y = 1). It can be seen that indeed, when B_s is increased, the reduction of the mass fraction due to an increase in the Kn number (from Kn = 0.1 to Kn = 1) is compensated, resulting into a similar profile. This confirms that by increasing the slip factor, which amounts to say altering the material properties of the scaffold pores to properties that reduce the interaction with the drug molecule solutions, compensates the reduction of drug release in nanopores.

Figure 9.  Mass flux as a function of time at x = 1 for different values of Kn and B_s = 0 (left) and 10 (right).

DownLoad: Full-Size Img PowerPoint

5.   Conclusions

This work presents new developments in Extended Non-Equilibrium Thermodynamics, proposing a model for studying drug uptake through cell membranes and drug release out of scaffolds. The theory is based on elevating fluxes to the same independent status of the usual variables. In this particular case, the mass flux has the same status as the mass fraction or concentration of the medication solution. This means that the mass flux follows its own evolution next to that of the mass fraction, being connected on an equal level via a species balance. By taking into account the relative importance of the relaxation time with respect to the diffusion characteristic time, and of the mean free path with respect to the system characteristic size, non-local effects are included in this way.

As for the diffusion through a cell membrane, it appeared that for sizes and relaxation times of the order of magnitude of or smaller than the mean free path (Kn = O(1)) and diffusion characteristic time (M = O(1)), respectively, leads to a lower drug uptake. This is due to the non-local effect, mentioned before. Non-Fickian profiles result from this, where it can be seen that at high non-locality a large mass fraction gradient is present close to the membrane entrance, but a nearly horizontal one near the exit. The latter contributes to the reduction of drug uptake. Furthermore, it has also been shown that at the membrane exit, a mass fraction jump (within the membrane) occurs that is not due to the boundary condition, but also due the non-locality effect. This is an effect that has not only been observed for heat transport ^[25], where it is already an established phenomenon, but also for mass transport ^[26]. The effect of the partition coefficient on the mass fraction is also studied, since they are typically introduced to account for the typical mass fraction jump (outside the membrane) observed in membrane diffusion studies, though still presenting experimental uncertainty to its precise value ^[26,27,28]. No particular effect has been observed other than that a higher partition coefficient causes a higher mass fraction jump at the membrane boundaries. If the medication could be chosen in such a way that the partition coefficient is high at the entrance and low at the exit, the efficiency of the medication could be increased if a higher concentration would lead to a better healing. Of course, mass conservation implies that this would not influence the total quantity of the uptake, but rather the speed of uptake, which could influence the kinetics of healing.

As for the drug release from the scaffold, the non-locality effect also showed to have a reducing impact on the release of drugs from the scaffold. It appeared that such a reduction could be countered by choosing the material of the scaffold in such a way that (besides biocompatibility, which is not a complicated task on its own) the interaction between the scaffold pore walls and the medication solution is weakened. This results into a higher slip factor enhancing drug release.

Future developments of this work include experimental demonstration, using the proposed model in this work. As quantitative agreement between experiments and models are difficult to find in the literature, comparison with experimental data should first be performed by nanomaterials of which material properties are well known and documented. As example, we can mention the delivery of nucleic acids with DNA nanostructures as new means of drug delivery by nanocarriers ^[29,30]. In the case of the nanocarriers, which can also be seen as scaffolds, the diffusion of the carriers through the cell membranes could be described by the set of Eqs (10)–(12) and (14)–(15). Subsequently, the delivery of the drugs out of the carriers into the cells could be described as release from a scaffold by Eqs (20)–(22). The difficulty lies in determining the interaction the medication can have with the nanocarrier, such as the DNA nanostructures ^[29]. In other words, it is important to know the slip length or slipping coefficient in Eq (17). Enhancement of drug release does not depend only on the release of the drug out of the carrier, but also the solubility of the carrier itself through the cell membranes. It is reported that DNA nanostructures can overcome the usually large circulation times of carriers ^[29,30]. The aforementioned discussion implies a necessity of measuring some material properties in order to use the presented model in this paper. In some cases, information about material properties could be obtained by in vitro experiments with a controlled experimental setup in a similar way performed by ^[31] for high-precision regulation of cell behaviour. This brings the intermediate step between modelling and full experimental demonstration to be identified by bio-mimicking cell membranes or scaffolds. Finally, such information can complete an experimental demonstration of our model.

Acknowledgements

Financial support from BelSPo is acknowledged.

Conflict of interest

The author declares there is no conflict of interest.