Surging footprints of mathematical modeling for prediction of transdermal permeability.

In vivo skin permeation studies are considered gold standard but are difficult to perform and evaluate due to ethical issues and complexity of process involved. In recent past, a useful tool has been developed by combining the computational modeling and experimental data for expounding biological complexity. Modeling of percutaneous permeation studies provides an ethical and viable alternative to laboratory experimentation. Scientists are exploring complex models in magnificent details with advancement in computational power and technology. Mathematical models of skin permeability are highly relevant with respect to transdermal drug delivery, assessment of dermal exposure to industrial and environmental hazards as well as in developing fundamental understanding of biotransport processes. Present review focuses on various mathematical models developed till now for the transdermal drug delivery along with their applications.

Introduction

Skin is the largest organ of the human body having a very complex structure. Due to unique structural and physicochemical properties, it is very different from other biological and microporous membranes. It consists of multi-layers including epidermis (thin outermost layer), dermis (a thicker middle layer) and subcutaneous tissue layer i.e. hypodermis (innermost layer). The skin performs three main functions i.e. protection, regulation and sensation. The regulatory function of the skin attracts the interest of scientists for developing formulations for skin [1]. Transdermal permeation occurs through three pathways namely: the diffusion through the lipid lamellae; the transcellular diffusion through the keratinocytes and lipid lamellae; permeation through appendages, hair follicles and sweat glands. The drugs should have sufficient lipophilicity to partition into SC but also should have sufficient hydrophilicity to pass through the epidermis and eventually through the systemic circulation. For most of the drugs, the rate determining step for drug transport is transit across the SC [2]. A large number of experimental and theoretical investigations have been carried out on the skin permeability. The prediction of percutaneous permeation is attracting attention of researchers in Cosmeceutical and pharmaceutical industry. Hence, development of mathematical models of epidermal and dermal transport seemed to be essential for the optimization of percutaneous delivery of drugs and for evaluation of their toxicity [3]. Mathematical models of skin permeability are highly relevant to the fields of transdermal drug delivery and in developing fundamental understanding of biotransport processes. Modeling of percutaneous permeation provides an ethical and viable alternative to laboratory experimentation. In vivo skin permeation studies are considered gold standard, but are difficult to perform and evaluate due to ethical issues and complexity of the process involved [4]. In vitro measurement of skin permeation can be done simply by using diffusion cell. Although it is easy and viable, this method is time consuming. In light of the above factors, research in developing mathematical modeling for transdermal drug delivery is at a high pace these days.

Mathematical models are the collection of mathematical quantities, operations and relations together with their definitions and they must be realistic and practical. The mathematical model is based on the hypotheses that consider mathematical terms to concisely describe the quantitative relationships. Many models have been proposed till now; however, earlier mathematical modeling was not in that much progress as in present scenario. In addition to establishing the required mathematical framework to describe these models, efforts have also been made to determine the key parameters that are required for the use of these models. The first contribution to mathematical modeling was given by Takeru Higuchi; a pharmaceutical scientist who applied physical and chemical principles to the design of controlled release devices in 1961 [5]. He proposed an equation exhibiting a considerable initial excess of undissolved drug within an inert matrix with film geometry allowing for a surprisingly simple description of drug release from an ointment base. The importance of mathematical modeling was more clearly understood by the year 1979. Categorically, mathematical models can be divided into empirical and mechanistic models. However, detailed number of mathematical models developed for analyzing and predicting data of transdermal studies are summarized in Fig. 1.

Fig. 1

Classification of various mathematical governing transdermal permeation models.

Empirical models

Empirical models are based on experimental data. These models are not based on physical principles and also not on assumptions made with respect to relationship between different variables. Empirical models are computer based modeling developed by Meuring Beynon in early 1980s. The main software used in empirical modeling is TKEDEN.

Multiple linear regression models

Multiple linear regression models help in relating two or more independent variables and a dependent variable by fitting a linear equation of the data obtained from the observations. A multiple regression model for independent variables x having parameters α for calculation of the dependent variable y is given bywhere,

i = 1, 2, 3…, n observations

Multiple regression analysis can also be done by least-square analysis model. Potts and Guy predicted skin permeability by using this model [6]. They obtained permeability coefficient data for transport of a large group of compounds through the mammalian epidermis, which was analyzed by multiple linear regression model based upon permeant size and octanol/water partition coefficient. These analytical data helped in predicting the percutaneous flux of pharmacological and toxic compounds entirely on the basis of their physicochemical properties.

Multiple regression analysis is one of the old techniques and is being used by several investigators till date (Table 1). The prediction of dependable variable is possible by using multiple independent variables. The major advantage of multiple regression models is that non-optimal combinations of predictors can be avoided. These models allow the examination of more sophisticated research hypotheses than is possible using simple correlations and it links various correlations with ANOVA models. This is an exceptionally flexible method. The independent variables can be numeric or categorical, and interactions between variables can be incorporated; and polynomial terms can also be included in the model.

Table 1

Applications of empirical models.

Reference	Experiment performed	Inference
Multiple Linear Regression Model
Potts and Guy [6]	Permeability coefficient data for transport of a large group of compounds through mammalian epidermisAnalyzed by multiple linear regression model based upon permeant size and octanol/water partition coefficient	Predicted percutaneous flux of pharmacological and toxic compounds entirely on the basis of their physicochemical properties
Sartorelli et al. [7]	Percutaneous diffusion of 16 compounds, eight of which were polycyclic aromatic hydrocarbons, six organophosphorous insecticides and two phenoxycarboxylic herbicides, were tested in vitro using monkey skin.Log octanol/water partition coefficient values were correlated with experimentally determined values of the permeability constant and lag time.	Precise values of permeability and other physicochemical parameters on the basis of log octanol/water partition coefficient and water values were predicted experimentally, using the algorithm derived from the multiple linear regression equation.A good correlation between percutaneous absorption data and physicochemical properties of industrial chemicals was established.
Models Based on Artificial Neural Network
Takayama et al. [8]	Applied an ANN system to a design of a ketoprofen hydrogel containing O-ethylmenthol (MET) to evaluate the promoting effect of MET on the percutaneous absorption of ketoprofen from alcoholic hydrogels in rats in vitro and in vivoThe amount of ethanol and MET were chosen as causal factors.The rate of penetration, lag time and total irritation score were selected as response variables.A set of causal factors and response variables was used as tutorial data for ANN and fed into a computer.	Nonlinear relationships between the causal factors and the response variables were represented well with the response surface predicted by ANN.Ketoprofen hydrogel was optimized using generalized distance function method.
Kandimmala et al. [9]	Developed to optimize a suitable vehicle composition to deliver melatonin via transdermal route.Solvents like water, ethanol, propylene glycol, their binary and ternary mixtures were used to increase the flux of melatonin and to reduce the lag time.Special quartic model (response Surface Method) and ANN were employed as prediction tools.	Water : Ethanol : Propylene glycol in the ratio of 20:60:20 had showed the best permeability (12.75 µg/cm2/h) and lag time 5 h.The various responses (solubility, flux, and lag time) summarized by response surface method and ANN with respect to vehicle composition helped in studying the inter-relativity between the responses.
Takahara et al. [10]	Optimization of ketoprofen hydrogel containing MET was done using ANN system.12 types of hydrogels were prepared using varying amount of ethanol and ketoprofen and study was performed in vivo in rats.Rate of penetration, lag time and total irritation score were chosen as response variables.For tutorial data for ANN, set of causal factors and response variables was used and fed into a computer.Optimization of the ketoprofen hydrogel was done according to the general distance function method.	Nonlinear relationships between the causal factors and the release parameters were represented well with the response surface predicted by ANN.Experimental results of rate of penetration and total irritation score coincided well with the predictions.It was inferred that the multi-objective simultaneous optimization technique using ANN was useful in optimizing pharmaceutical formulas when pharmaceutical responses were nonlinearly related to the formulae and process variables.
Obata et al. [11]	The effect of 35 newly synthesized MET derivatives on percutaneous absorption of ketoprofen was investigated in rats by using ANN system which helped in understanding the relationship between the structure of compounds and promoting activity i.e. structure–activity relationship.An enhancement factor equal to the ratio of rate of penetration with enhancer to the rate of penetration without enhancer and total irritation score were chosen as response variables.Factors like log P, molecular weight, steric energy, van der Waals area, van der Waals volume, dipole moment, highest occupied molecular orbital and lowest unoccupied molecular orbital were used to find out the structural nature of cyclohexanol derivatives.	The experimental values of enhancement factor and total irritation score get coincided well with the predicted values by ANN design.As the contribution of log P in predicting the enhancement factor was almost equal to 50%, it has been inferred that promoting activity of these compounds is mostly affected by their lypophilicity only as compared to other physicochemical properties.
Degim et al. [12]	ANN analysis to predict the skin permeability of 40 xenobioticsPermeability coefficient was the key parameter.Used a previously reported equation for prediction of skin permeability by using the partial charges of the penetrants, their molecular weight and octanol water partition coefficients	Correlated experimental and predicted permeable coefficient values from literature and the regression value was found to be 0.997.Advantage-ANN model developed does not require any experimental parameters; it potentially provided a useful and precise prediction of skin penetration for new chemical entities in terms of both therapy and toxicity.Reduced the need of performing penetration experiments using biological or other model membranes

Besides these advantages, the model is associated with certain limitations like unstable regression weights and poor repeatability etc. In addition, large samples are desirable and, although no exact guidelines are specified regarding this, however, number of samples should considerably exceed the number of variables. Another limitation of MLR is its sensitivity to outliers like if most of our data is in the range of “20, 50” on the x-axis, but we have one or two points out at x = 150, this could significantly swing our regression results. One more issue is overfitting. It is easy to overfit the model such that the regression begins to model the random error in the data, rather than just the relationship between the variables. This most commonly arises when too many parameters are compared to the number of samples. Linear regressions are meant to describe linear relationships between variables. So, if there is a nonlinear relationship, then it is not a best fit model. However, in certain cases it can be compensated by transforming some of the parameters with a log, square root, etc.

Models based on artificial neural network (ANN) system

ANN is a computer-based system and inspired from a simple neural structure of the brain where a number of neurons/units are interconnected in a net-like structure. It is also known as feed-forward layered neural network (FFLNN) consisting of an input layer, an output layer and many intermediate layers or hidden layers. Each unit in a layer is influenced by other units in adjacent layers. Their values of connections or weights affect the degree of influence of neurons. Till now, this approach has been successfully used to predict the octanol/water partition coefficient (ko/w), oral bioavailability (BA) of drugs for analysis of clinical pharmacokinetic data and ultimately in designing of pharmaceutical formulations [13]. ANN systems are of many types varying from one or two layered single directional to complicated multi-input directional feedback loop layers. A highly trained person is required for operation and calculating such models.

ANN models can be classified into three variables such as various significant formulations and categories based on their functions, associating networks, feature extracting networks and non-adaptive networks. Associating networks are employed for data classification and prediction of need input (independent variable) and correlated output (dependent variable) values to perform supervised learning. Feature-extracting networks, which are used for data (ANN) dimension reduction and need only input values to perform unsupervised or competitive learning. Non-adaptive network needs input values to learn the pattern of the inputs and reconstruct them when the computer is presented with an incomplete data set [14].

ANN models have certain limitations also. The ANNs modeling is an alternative to conventional modeling techniques. ANNs have gained application in modeling the process which cannot be tackled by classical methods. The ANNs do not require special software or computer as they can be described using simple function of computers. These systems are better than mathematical models, e.g. response surface methodology has the potential to solve and recognize problems involving complex patterns. ANNs can predict chemical properties of compounds in a better way than MLR, e.g. solubility of APIs. They cannot be used to elucidate the mechanistic nature of the correlation established between the variables. A formulator may need a lot of proficiency to obtain a reliable ANN model. The front end of the work such as experimental design and data collection may be more time consuming than the traditional approach used by experienced formulation scientists [7].

The utility of ANN models lies in the fact that they can be used to infer a function from observations, especially by hand in those cases where the complexity of the data or task makes the design of such a function impracticable by hand. ANN can be used in regression analysis, fitness approximation, data processing, robotics, prosthesis and can also be used in computational neuroscience. Table 1 summarizes the work done by different scientists on empirical models.

Membrane transport models governing iontophoretic delivery

The designing of information for any model is an attempt to describe the complexity of a real biological system. It also tests that the particular model deals with the complexity of the model with how much accuracy. The research involving biological systems is associated with increased level of complexity. This complexity arises due to the anatomy of the system, feedback-control loops, biochemical reactions within the tissue and due to the presence of specialized transport mechanisms for specific substrates. The complexity present in the system often supersedes the prediction of its behavior. Hence, artificial aids have been used to understand such complexities and this is done by developing different mathematical models for membrane transport. A number of models have been developed and used to predict the transdermal permeation profile after application of iontophoresis. Iontophoresis is a technique that involves the application of a small voltage across the skin to drive ions into and across the membrane [15] (Fig. 2). Since ion flow is facilitated by transportation through aqueous pathways, the physicochemical properties that decrease the passive diffusion through the intercellular lipidic space favor electrically-assisted delivery. The Iontophoretic transport rate depends on the intensity, duration and profile of current applied. Iontophoresis can be used for the controlled delivery of therapeutic molecules including peptides and proteins. Researchers have compared it with a “needle-less” infusion pump. A number of factors that significantly affect iontophoretic drug delivery are ion composition, solute size, charge, solute mobility, total current applied, presence of extraneous ions, epidermal permeation selectivity, pore size, etc. No single model can integrate all these determinants; hence a single model cannot satisfy all the practical conditions. However, models have been designed in order to integrate maximum of these determinants. Still, the work is going on in this area of research.

Fig. 2

Ionic movement across the membrane during iontophoresis.

In Table 2, some of the models related to iontophoretic delivery are summarized along with their applications.

Table 2

Applications of membrane transport models governing iontophoretic delivery.

Reference	Experiment performed	Inference
Poisson–Nernst–Planck (PNP) Model
Shrinivasan et al. [16]	Transdermal flux of triethlyammonium (TEA) across hairless mouse skin (HMS) was investigated after the voltage drop and the water flux for model charged solutes at steady-state, using the four-electrode potentiostat system.The method of correcting for the skin damage effect was introduced.The synergism in iontophoresis and pretreatment enhancer (ethanol) was also investigated for delivering a high molecular weight polypeptide.	The contribution of water transport on solute flux was observed to be lower than the contribution due to the applied voltage drop.The theoretical predictions of iontophoretic flux were higher than the experimental observations.Pretreatment with ethanol in combination with iontophoresis influences the permeability coefficient of insulin.
Kasting and Bowman [17]	Human allograft skin was immersed in saline buffer and direct current–voltage relationships and sodium ion transport measurements were determined using diffusion cell and four terminal potentiometric method.	Sodium ion permeability coefficients by this method were less as compared to permeability coefficients of sodium ions in human skin in vivo.Current–voltage relationship in the tissues was found to be time dependent and highly non-linear.Resistance of skin decreased with increase in current or voltageFlux was found to be ~3–5 folds more than the predicted values.
Kasting et al. [18]	The validity of Nernst–Planck equation was tested for homogenous membrane under the constant-field for steady- state and unsteady state through skin.Validity was done by observing the iontophoretic transport of a negatively charged bone resorption agent, etidronate disodium (ethanehydroxydiphosphonate, EHDP) across excised human skin at different voltages and currents.	The iontophoretic transport results were highly variable under constant voltage or constant current which invades the passive skin transport area.
Nernst-Planck (NP) models with convective flow
Tojo [19]	An iontophoretic model based on time-dependent drug binding and metabolism as well as the convective flow of solvent was developed.Solutions were obtained by numerical integration of the resulting partial differential equation under the constant field approximation.Effects of mode of application, electric potential, diffusion coefficient of the drug, skin-drug binding and convective flow across the skin caused by the electric field on skin permeation and plasma concentration were also produced, and these data were also compared with in vitro transdermal iontophoretic data of a polypeptide.	Mode of application affects the permeation profile and plasma concentration profile quantitatively.Iontophoretic transdermal delivery is effective for large molecules such as peptides and proteins that penetrate into skin with great difficulty by conventional passive diffusion.
Shrinivasan and Higuchi [20]	Effects of the applied electric field and convective solvent flow on the permeant flux were investigated.Predictions of this model were compared with a model in which no convective solvent flow term was considered	Provided a better detailed framework to decouple and understand the interactions of the applied field and the solvent flow effectsThe effect of the convective solvent flow on the iontophoretic flux was found to be inversely related to the molecular size (diffusion coefficient) of the permeant.
Hoogstraate et al. [21]	Studied the iontophoretic increase in transdermal transport of leuprolide in vitro.An exhaustive investigation was done to discover the mechanisms of the inconsistent behavior of the positively charged peptide.Used a model membrane as well as human skin	Adsorption of leuprolide on to the negatively charged membrane leads to a change in the net membrane charge and therefore changing the direction of the electroosmotic flow.Due to reversal in the direction of the electroosmotic flow, the convective solvent flow has been hindered rather than assisting the flux of a positively charged permeant and assisted the iontophoretic flux of negatively charged permeants.
Imanidis and Leutolf [22]	Analyzed the experimental increase in flux of an amphoteric weak electrolyte measured in vitro using human cadaver epidermis at a voltage of 250 mV at different pH values.	The shift of pH in the epidermis when compared to the bulk was caused by the electrical double layer at the lipid-aqueous domain interface which was evaluated using the Poisson–Boltzmann equation.The enhanced flux depends upon factors such as applied voltage, convective flow velocity due to electroosmosis, ratio of lipid to aqueous pathway passive permeability, and weighted average net ionic valencies of the permeant in the aqueous epidermis domain.The model can provide a good quantitative insight into the reciprocation between different phenomena and permeant properties influencing iontophoresis.
Ferreira et al. [23]	Presented a multi-layer mathematical model using NP equation with convection–diffusion process to describe the transdermal drug release from an iontophoretic system.The stability of the mathematical problem is discussed in two scenarios i.e. imperfect and perfect contact between the reservoir and the target tissue.An accurate finite-difference method is proposed to explain the drug dynamics in vehicle and skin layers coupled during and after the electric administration.	This multi-layer model was developed to clarify the role of the applied voltage, the diffusion of drug, the conductivity of the skin, and the systemic absorption.This model was found to be a simple and useful tool in finding new delivery strategies that ensures the optimum and localized release of drugs for a prolonged period of time.
Roberts et al. [24]	Developed an integrated ionic mobility-pore model for iontophoresis of epidermis using two types of models which are free volume type model and pore restriction type modelThis model was developed for finding out some parameters of iontophoretic model like the solute ionic mobility in the aqueous solution present in the pore and in the donor solution, the effect of pore size restriction on iontophoretic solute and the effect of ionization of partial solute.Used a number of solutes and developed a model in which examination of the determinants of transport for individual solutes was done where ionic mobility and size are kept as constant parameters	It has been found that iontophoretic transport of a number of ionizable solutes depends on pH and the extent to which solute interacts with the pore wall depends on the fraction of ionization.

Poisson-Nernst–Planck (PNP) model

The molecular mechanism of ionic movements through transmembrane channels is one of the most interesting studies in biophysics, which is of great importance in living cells. Now the question arises what ion channels are. Ion channels are pore-forming proteins found in cell membranes which allow only specific ions to pass across the membranes and help in maintaining optimum ionic composition. Hence, they help in regulating cellular activity via maintaining the ionic flow [25] and acts as vital elements in maintaining many biological processes like excitation, signaling, gene regulation, secretion and absorption [26]. Therefore, ion channels are very important for cell survival and function. A number of theoretical and computational approaches have been developed over the past few decades to understand the physiological functions of the ion channels, and PNP model is one of the most popular approaches among all.

The PNP model is based on a mean-field approximation of ionic interactions and continuing description of electrostatic potential and concentration. This model provides qualitative descriptions and quantitative predictions of experimental observations for the ion-transport models in many fields like nanofluidic systems and biological systems. In this theory, the Poisson equation is used to describe the electric field in terms of electrostatic potential. The electrodiffusion of ions in terms of ion concentration is described by Nernst–Planck equation. Planck has given a general mathematical solution to problems related with electrodiffusion of ions [27]. Ion mobilities present in a solution were related to their diffusion coefficients by Nernst [28]. Using these results, Planck has considered particular solutions to the following set of differential equations which are now known as Nernst–Planck flux equations:where,

j = Flux of ionic species i in a convection-free fluid or membrane

D = Diffusion coefficient

c = Concentration of ionic species present in the fluid or membrane

z = Charge present on ionic species

F = Faraday constant

R = Gas constant,

T = Absolute temperature

= Electric potential.

This equation describes one-dimensional transport in dilute solutions [29], [30]. Any solution which describes ion transport in membrane must satisfy the Poisson equation in one-dimension if it is believed to be harmonious with electrostatic theory:where,

∈ = Permittivity of the membrane

is the local space charge density (e is the electronic charge).

Equation (2) describes the statement of conservation of charge for a continuum.

Equation (1) and Equation (2) specify the ionic transport within a system in which solvent velocity is zero (i.e. unstirred system). Both of the above mentioned equations are nonlinear differential equations which must be solved numerically. In certain situations, simplified approximations can be made and analytical solutions to the problems can be generated using these equations e.g. in case of very thin or very thick membranes or equal ion concentrations on both sides of the membrane.

PNP model is a very successful model as it is a good predictor of ion-transport phenomenon in biological channels for non-equational systems. However, it has some limitations also, e.g. it neglects the finite volume effect of ion particles and correlation effects and is of much importance with respect to ion transport in confined channels [31]. Another major drawback of PNP model is its high computational cost. Computational cost increases with increase in number of ion species. The number of equations, number of diffusion coefficient profiles to be determined depend on the number of ionic species in system as each ionic species corresponds to one Nernst–Planck equation and to one diffusion coefficient profile. Hence, a complex system with multiple ion species will cost very high. For such systems, Nernst–Planck equations are substituted with Boltzmann distributions of ion-concentrations.

PNP theory cannot be used to explain the transport in channels as average number of ions in the channels is comparable to the fluctuations in size and hence the concept of concentration gradient does not hold true. However, PNP theory explains the saturation of ionic flux as a function of ionic concentration in the solution adjacent to the biomembranes so as to attain the fixed membrane potential.

Poisson–Boltzmann–Nernst–Planck (PBNP) model

To solve all the problems in solving multiple ion species in a complex system, an alternative model was suggested and his model is popularly known as PBNP model. This model is derived from total energy functions by using the variational principle. By simply modifying the PNP model like including the stearic effect of ionic transport, a qualitative result has been observed in PNP computational results [32], [33].

After the development of univalent ions by Moore [34], dimensionless variables ζ = x/h, (or , where k is Boltzmann's constant), , where refer to the negative ions, c refer to the positive ions and is the average concentration of total ions in the membrane. Now the Poisson equation can be written as:where , the Debye length, is the characteristic width of the space charge layers at the boundaries of the membrane. The quantity p − n shows the excess of positive charge over negative charge at any point. Different solutions can be generated depending upon the width of the membrane, h and λ.

Despite of the great success of PNP model, it has some limitations also, such as neglecting the finite volume effect of ion particles and co-relation effects. Such limitations become very important in highly confined ion channels [35].

Nernst–Planck (NP) models with convective flow

While considering the formal description of steady-state diffusion of charged permeants under an electrochemical potential gradient, Nernst–Planck equation is most commonly used as its starting point [36]. The Nernst–Planck equation is an example of a single flux model, where the conjugate driving force (the electrochemical gradient) is assumed to drive the flux completely and the non-conjugated flows or forces are neglected. In short, in Nernst–Planck model, the transport of a permeant due to convective flow of solvent is assumed to be negligible. However, an electrically driven flow of ions across a membrane carrying a net charge induces a convective flow of solvent through a phenomenon known as electro-osmosis [37]. In case of charged membranes, the overall electroneutrality requires an excess of mobile ions of opposite charge within the membrane. In the case where cations and anions present in the system have comparable mobilities, the species with higher concentration will carry more current when an electrical field is applied on the membrane. This can lead to a net flow of solvent in the direction of the carrier having higher concentration. The effect of convective solvent flow can be shown by adding to the right hand side of the Equation (1). Then the NP equation for convective flow will become:where,

v = Solvent velocity

c = Solute concentration and

σ = Reflection coefficient at the membrane-solution boundary

Kinetic models of ion transport

Iontophoresis has received a great attention in recent years for delivery of peptides and other poorly bioavailable drugs through the skin [38]. Many theoretical approaches have been reviewed toward the quantitative description of this phenomenon [39], [40]. As we have discussed earlier, most of the models are based on Nernst–Poisson equation obtained under constant field approximation or extended to include electro-osmotic effects for convective flow of solvents by including solute-solvent coupling terms. Also, in case of thick membranes surrounded by solutions of different ionic strengths, the electroneutrality approach (Planck's approximation) has been considered best. The problem with using all these models is that none of them relate the exponentially non-linear and slightly asymmetric current-voltage observed in direct current experiments with excised skin. Hence, kinetic models have been developed for finding solutions to such problems. As we know, Nernst–Planck equation models treat the membrane as a continuum. The conductance of the membrane is due to the space charge layers and mobilities of ions within the membrane. Boundary conditions must be applied to relate the ionic concentration in the membrane as well as that of outside the membrane. The electrochemical potentials of each ionic species on opposite sides of the membrane–solution interface can be related with equilibrium conditions using following equations:

However, this equilibrium may not be satisfactory in the case where current flows through a membrane after imposition of electrical potential. In such cases, kinetic model can be used to provide the boundary conditions for ion transport. Hence, the equation generated will be in the form:where and are rate constants given by:

In case of static barriers, i.e. independent of concentration and voltage; the transport properties are almost similar to that of equilibrium system described in Equation (1), (2) and (5). However, a great difference in transport properties is observed if the magnitude of the potential barrier is not static and is dependent on applied voltage and concentration. This was first noted by Burnett [39]. A model for ion transport in lipid bilayer membranes with respect to overpotential analogous was described by Bender [41] as the following equation:where z is charge of electrolyte in which membrane is dipped having concentration c and c in opposite sides. Further, Butler-Volmer theory of electrode reaction kinetics [42] may be applied to ion transport in membrane, where system is dealing with diffusion and kinetically limited rates of charge transfer:where,

c(0,t) and c(0,t) are time-dependent functions denoting the electrolyte concentration at the surface of the membrane; α denotes the symmetry measure of the energy barrier (normally, α = 0.5 for symmetrical barriers).

Equation (10) is used where a depletion layer develops at the solution–membrane interface due to finite diffusiveness of the ions. However, Equation (11) applies when this process is not taken into consideration. This concept of voltage-dependent ion transport has been applied to many membrane transport problems [41], [43].

The electrochemical processes at the interface of electrolyte solution in water and an organic liquid utilized the experimental methods of double layer studies and ion transfer kinetics. The study elucidated that equilibrium potential, double layer structure and kinetics of ion transfer are dependent on each other [44].

Hindered transport model

In order to provide a rationale account for the complex structure of the skin, many scientists have used Hindered Transport Theory (HTT) to develop mathematical modeling for transdermal iontophoresis [45]. Hindered transport theory is used as an extension of the Nernst–Planck equation in the studies where the effects of constrained flow geometries and electrostatic interactions on the fluid are considered. In this model, the hypothesized aqueous pores of SC are considered as integral structures of the skin through which the transport occurs. These pores can be temporary channels also that are induced after the application of current.

The hindered transport model was initially developed to characterize flow through long, narrow passages, such as capillaries and straight channeled porous membranes [46], [47], [48], [49]. The overall resistance to flow through the membrane was strongly affected by interparticle interactions and particle–wall interactions. The same types of interactions are thought to be present in large channels; however, in their case the relative influence on the hydrodynamic flow profile could also be neglected.

The flux in HTT can be represented by Equation (12)where,

ε = Tortuosity factor which is lumped together with void fraction

H = Hindrance factor for diffusion and migration

W = Hindrance factor for convection.

On a molecule scale, H accounts for steric and long range electrostatic interactions and W accounts for the enhanced hydrodynamic drag on particles caused by the presence of the pore wall. According to Anderson and Quinn [46], the hindrance factors, i.e. H and W for spherical particles are given below:where the independent variable λ is defined by Equation (16)where is the particle radius and is the pore radius.

Refined hindered transport model governing ionic mobility with respect to pore wall

A major benefit of the hindered transport model is that the negative charge of the skin could be considered in this model. An assumption is generally considered that the charge is present on the surface of pores which are cylindrical in shape. Due to the negative charge on the walls of the pores, positive charge get develops in the electrolytic solution adjoining to the pore surface which can get diffused through these pores. In this way, electroneutrality is maintained in the skin.

The thickness of this diffusion region can be estimated by the Debye screening length. The equation for calculating Debye length is:where,

ε = Solution permittivity, or dielectric constant

R = Universal gas constant

T = Absolute temperature

F = Faraday's constant

z = Bulk solution concentration

= Number of charged ions in bulk solution

A force will be exerted on the volume of ionic solutions present in this diffusional region when an electrical field is applied [50]. If this electric field is perpendicular to the pore walls, bulk fluid flow starts. This phenomenon of bulk fluid flow when an electric field is applied is known as electroosmosis.

As skin is having a complex morphological structure, it seems to be very unrealistic that transport of charge occurs through straight channeled pores. A refined hindered transport model of transdermal iontophoresis which include solute interactions with pore walls was suggested [24]. In this model, they considered the effects of partially ionized solutes and irregularly shaped particles. This model has developed an expression which related the molecular volume of a drug with iontophoretic rate. In this refined model, flux of a species is calculated as:where,

= Overall iontophoretic permeability coefficient

c = Solute concentration

It was assumed that the transport of charged molecules across the human skin could be increased by applying high-strength electric field. This was theoretically observed by Edward and his co-workers [51] in terms of electroporation of lipid bilayers present in SC above the voltage of transbilayer for which electropores have been observed in single bilayer membranes. Electroporation includes the formation of ephemeral aqueous pathways in lipid bilayers of a brief electric pulse. This phenomenon occurs when a voltage reaches 0.5–1 V for short pulses across a lipid bilayer. The transdermal molecular flux was predicted at two electric field conditions, according to the size, shape and charge of the transporting molecule. These two electric field conditions were either a small electric field, i.e. transdermal voltage <<100 V or sufficiently large electric field i.e. transdermal voltage >100 V. At small field strength, charged molecules are transported through shunt routes of the skin and at large strength electric fields, charged molecules accessed a transcorneocyte pathway and transbilayer transport occurs through electropores of lipid bilayers. An increase in transport was observed during skin electroporation.

The limitation of this model is that the estimation of flux of preferred drug cannot be done because the relationship between drug physicochemical properties of drug and pore size has to be determined experimentally.

Nonequilibrium thermodynamic models

Numerous mathematical models have been developed for molecular transport across biological membranes based on nonequilibrium thermodynamics [52], [53]. Nonequilibrium thermodynamic models of transdermal iontophoresis are fascinating because upon application of current model, skin will not be under equilibrium conditions. Hence, only a limited number of attempts have been made to apply the infrastructure to transdermal iontophoresis [54], [55], [56]. The application of nonequilibrium thermodynamics for modeling biological membrane transport was spearheaded by Kedem and Katchalsky [57], [58], [59]. The specific interactions between the membrane and the electrolyte solution components were being accounted in the investigation. In this work, the solvent has been considered as a diffusing species and the undeniably bulk fluid flow was also considered as is in case of biological membrane transport [57]. The data treatment represented an abandonment of the classical Nernst–Planck formalism governing diffusion processes.

Nonequilibrium or irreversible thermodynamics are based on the assumption that there is a direct relation between the forces and fluxes of a given system. This assumption also holds for equilibrium thermodynamics, however, the approach for defining the flux equations is different. For example, the electrochemical potential gradient is considered as the driving force for the flux in the Nernst–Planck formalism. However, in nonequilibrium thermodynamics, the forces and fluxes are constrained by the dissipation function. The integrated form of the dissipation function anticipation function is defined according to Equation (19);where,

= Integrated form of the dissipation function

= Rate of entropy generation per unit volume of the membrane

T = Absolute temperature, which is placed outside the integral because biological systems are essentially isothermal.

The integrated form of the dissipation function is usually applied for membrane transport because it is difficult to assess the local forces within the interior of the membrane.

The fluxes are related to the dissipation function according to Equation (20);where,

Ji's = Fluxes

Xi's = Driving forces for the fluxes.

Equation (20) states that the rate at which the fluxes contribute to the entropy of the system is proportional to the driving force. The general form of species flux for irreversible thermodynamics, subject to the constraint defined by Equation (20), is described bywhere , is the phenomenological coefficient. The phenomenological coefficients are proportionality constants which represent the contribution to the species fluxes from a given force.

The Q is included to account for the interaction of species i, with all other components in the solution. The approach accounts for interactions between the solvent and the various solute molecules.

Macroscopic modeling

Percutaneous permeation occurs on a spectrum of pecking order or scale. Different phases of penetration have distinctive length scale. The quantitative consideration of transport phenomena is making a broad entry into the modeling of cellular processes at several scales as shown in Table 3.

Table 3

Different length scales representing skin as a barrier.

Level	Scale	Possible area/expectational level	Length
Level L1	Macromolecular level	Lipid bilayer	1–10 nm
Level L2	Cellular and sub-cellular level	Reference cell (periodic)	1–10 µm
Level L3	Membrane level	Penetration amount depending on depth of cell membrane in diffusion cell experiment	0.1–1 mm
Level L4	Compartmental level	Amount penetrated in different cell experiments per compartment or in body	>1 cm

Models developed at the membrane (L3) and compartmental (L4) levels are termed as macroscopic models. The macroscopic scale is defined as a mass or layer of many cells that appears as an effective homogeneous continuum. Traditional pharmacokinetic models at this scale have often represented such cellular masses in terms of well-mixed compartments that absorb or generate bioactive molecules and exchange them with other compartments. These models are based on laws related to rate and are expressed in terms of one or more phenomenological biochemical rate coefficients. These models relate the percutaneous permeation of skin and the vehicle as a consequence of a chronicle order of penetrant i.e. penetrant from vehicle to SC or deeper skin layers. Here, diffusion occurs either by a single pathway or as a result of a combination of multiple independent pathways. Macroscopic models represent multilamellar structure of skin where each layer exhibits diffusion coefficient and partition coefficient as its distinctive properties that are ultimately used to quantify solution of a solute.

Macroscopic transport is generally an outcome of an intricate interplay between physical and chemical processes on the microscopic scale at which the individual cells and cellular membranes are detectable as the microscopic elements are also comprised of rich substructures which determines their transport properties.

The macroscopic models are classified as a one-dimensional model (concentration in model membrane is related to the penetration depth) and two-dimensional model (considering drug transport through a membrane with respect to more than one membrane).

One-dimensional models

Three types of one dimensional models are described below and their applications are described in Table 4.

Table 4

Applications of one-dimensional macroscopic.

Reference	Experimental work	Inference
Models based on diffusion
Anissimov and Roberts [60]	A diffusion model was developed to evaluate the percutaneous absorption of a solute as a mean of constant donor concentration with a finite removal rate from the receptor due to either perfusion or sampling.The resistance of viable epidermis and resistance of donor SC interface were considered.Simulations of solute flux and cumulative amount absorbed and percutaneous absorption were calculated using numerical inversions of Laplace domain solutions.	The percutaneous studies were affected by experimental protocol, in-vivo perfusion conditions, perfusate flow-rate, finite receptor volume, resistance of viable epidermis/aqueous diffusion layer on percutaneous absorption kinetics and dermal concentrations.Limitation of this study was the inability to derive a simple Laplace equation.
Anissimov and Roberts [61]	Developed the model for investigating the effect of finite amount of solute on percutaneous absorptionNumerical invasion of Laplace domain solutions was used for data modeling.	Flux was found to be dependent on sampling rate of the drug.
Anissimov and Roberts [62]	Developed an another diffusion model for percutaneous permeation/penetration and desorption to find out the effects of varying diffusion and partition coefficients in stratum corneum over the diffusional path lengthUsed Laplace domain solutions to calculate steady-state flux, lag time and mean desorption time	Desorption model was found to be dependent on heterogeneity of diffusion and partition coefficients.
Kasting et al. [63]	Developed the model for taking into account the position of an arbitrary dose of a (potentially) volatile compound meant to be applied to the skinVehicle dynamics was neglected.Model was restricted to diffusion of a single volatile component only.	In the cases where dose was less than that required to saturate the upper layers of the SC, the shape of the absorption and evaporation profiles were found to be independent of the dose.If the dose was greater than that required to saturate the upper layers of the SC, the absorption and evaporation approach steady-state values with increase in dose.
Kruse et al. [64]	Combination of finite and infinite dose data measured of four compounds having different lipophilicities was utilized to obtain a new dermal penetration data.Two one-dimensional diffusion models were used.	Reproducible parameters of the absorption process were given successfully for finite and infinite dose by both types of models.Experimental permeability of a more lipophilic compound was found to be lesser than the theoretical permeability.Permeability by this model could provide more accurate and realistic values which can be used in QSARs and other applications.
Anissimov and Roberts [65]	Developed another diffusion model to show the effect of a slow equilibration/binding process within SC on absorption and desorption from SCDiffusion model solutions were used to derive the steady-state flux, lag time and mean desorption time for water in SC.	The effect of slow equilibration was less on the amount of solute absorbed than the amount of solute desorbed.
Models based on adsorption
Chandrasekaran et al. [66]	The sorption and rate of permeation of scopolamine in human skin was measured as a function of drug concentration in aqueous solution contacting the SC.Sorption occurred by both ordinary dissolution and binding of penetrant to immobile sites in the membrane with which the experimental sorption isotherm could be concluded, and also the disparity between steady state and time lag diffusivities could be resolved.	Sorption isotherm was found to be nonlinear.The apparent penetrant diffusivity computed from steady state permeation data was found to be greater than that estimated from unsteady state (time lag) measurements.It was concluded that sorption process is a simple model which conjured the coexistence of dissolved and mobile sorbed molecules in equilibrium with site bound and immobile molecules within the membrane, quite accurately correlated experimental sorption data and transient transport measurements.
Kubota and Twizell [67]	Examined the nonlinear mathematical model developed by Chandrasekaran and coworkers to oversee the pharmacokinetic profiles in percutaneous drug absorption (Chandrasekaran et al., 1980)A nonlinear partial differential equation of parabolic type was generated and developed a number of finite-difference methods for the numerical solution of the associated initial/boundary-value problem.Numerical differentiation and integration procedures were combined to oversee the cumulative amount of drug eliminated into the receptor cell per unit area as time increased.	The use of the equation for the simple membrane model to estimate the permeability coefficient and lag time is mandatory even if the system was described by the dual-sorption model, provided cumulative amount versus time data collected for a sufficient long time are used.Lag time was dose-dependent and decreased with increasing donor cell concentration.Permeability coefficient in the dual-sorption model remained constant and was found to be independent of the concentration of the donor cell.
Kubota et al. [68]	The nonlinear in vitro percutaneous permeation kinetics of Timolol was investigated with human cadaver skin.Several aqueous Timolol concentrations were kept in the donor cell for 25 h during the permeation studies.Analyzed the lag-time changes with a newly proposed method with the lag-time prolongation factor which was calculated from the parameter values obtained in the sorption isotherm studyDiffusion parameter for the mobile solute was determined by fitting the data of the cumulative amount excreted in the receptor cell versus time to the numerical solution of the dual sorption model.	Dual sorption model predicted that the plots of the amount of Timolol per unit area of epidermis versus aqueous Timolol concentration in equilibrium were curvilinear in the sorption isotherm (equilibrium) study.The same magnitude of the ratio of the lag-time to lag-time prolongation factor showed that the dual sorption model explains the nonlinear percutaneous permeation kinetics of Timolol.The numerical data of diffusion parameter was compatible with predicted values in the dual sorption model.
Gumel et al. [69]	Pade approximant (approximants derived by expanding a function as a ratio of two power series and determining both the numerator and denominator coefficients) (Gumel et al., 1993) had been used to develop the numerical method for the exponential term.The solution vector had been obtained via a two-stage sequential process by factorizing the rational approximant into its linear factors.	The non-linear, second-order type of parabolic partial differential equation was transformed in dual-sorption model.
Models based on metabolic activity in the living epidermis and dermis
Yu et al. [70]	A physical model approach was developed for the topical delivery of a vidarabine ester prodrug.Included modeling, theoretical simulations, experimental method development for factoring and quantifying parametersEmployed the deduced parameters to determine the steady-state species fluxes and concentration profiles in the target tissueHomogeneous enzyme distributions and constant diffusivities in the membrane were assumed for the simultaneous transport and bioconversion of the topically delivered prodrug.	Results provided the prevailing levels of the prodrug, the drug, and the metabolite at the target site and the transport rates of all species into the bloodstream.
Sugibayashi et al. [71]	The simultaneous skin transport and metabolism of ethyl nicotinate was analyzed.Many permeation studies were done to ethyl nicotinate or its metabolite nicotinic acid with or without using esterase inhibitor.Partition coefficient of ethyl nicotinate from the donor solution to the SC and diffusion coefficients of ethyl nicotinate and nicotinic acid through SC and the viable epidermis and dermis were determined.Michaelis constant and maximum metabolism rate were calculated from production rate of nicotinic acid from different concentrations of ethyl nicotinate in the skin homogenates.Fick's second law of diffusion and law of Michelis-Menton metabolism were used to analyze the data.	Steady-state fluxes of ethyl nicotinate and nicotinic acid were found at all the concentrations of ethyl nicotinate.Theoretically, the steady-state fluxes were calculated using Fick's second law of diffusion.Experimental data were in close proximity to the theoretical data.
Boderke et al. [72]	Developed a physical model to relate Fickian diffusion and concurrent Michaelis-Menten metabolism of drugs in the viable epidermis of human skinSubstrate concentration profiles were generated numerically within the metabolizing tissue and the resulting donor-to-receiver substrate fluxes through the tissue for various mass transport and metabolism parameters.Parameter required for those calculations were extracted from independent experiments.	Numerical predictions and experimental data were in excellent agreement.Principal validity of the physical model was proved by correlating independently validated parameters and numerical fits.Numerical simulations and theoretical derivations illustrated the kinetic impact of the factors involved i.e. the diffusion coefficient, substrate donor concentration, substrate partition coefficient, tissue thickness, and maximum metabolic rate on drug permeation.Advantage of this physical model is that the therapeutic or toxic substrate and metabolite levels in tissues may be calculated and the potential of a substrate to pass a metabolic barrier may be predicted.Since metabolism is a general feature of viable tissue, this model may be useful for the estimation of permeation and concurrent metabolism in other tissues and for other substrates as well.
Yamaguchi et al. [73]	Method was developed to analyze the simultaneous skin permeation and metabolism of 22-oxacalcitriol (follows a complicated topical metabolic pathway) by observing skin permeation of only unchanged 22-oxacalcitriol through excised rat skin.In vitro permeation experiments carried out using Franz-type diffusion cells.Time courses of unchanged 22-oxacalcitriol amount in ointment, skin, and receptor fluid were observed and these data were fitted to diffusion equations to obtain permeation parameters and the metabolic rate.	Fitting curves of the skin permeation profile obtained by the model were sufficiently close to observed data of unchanged drug amounts in ointment, skin, and receptor fluid. Due to this, calculation of skin permeation parameters i.e. partitioning, diffusivity, and metabolic rate of the drug was possible without requiring metabolic information.

Models based on diffusion

A prerequisite requirement in dosage form designing is the ability to quantitatively study the diffusion and membrane metabolism of drugs. Unfortunately, only few simple analytical solutions are available to model these processes. The predictive models to quantitate drug permeations were initiated by Scheuplein and workers [74], [75], [76]. They considered the skin barrier in terms of the fractional area and absolute thickness and at the same time diffusivity and partition coefficient were recognized as the important flux-determining properties.

The relatively simple mathematical equations can be used to quantitatively describe drug release model. The system which comprises diffusion controlled delivery can be of great use in speeding up the product development if applicable to a specific type of controlled drug delivery systems as it may allow in silico simulations of the effects of formulation and processing parameters on the resulting drug release kinetics. Additionally, more understanding of the underlying drug release mechanisms could be obtained from this model. However, caution should be taken so that none of the assumptions on which equations are based are violated.

Models based on adsorption

Generally, Dual Sorption Theory (DST) has been used to describe the process of permeation with respect to adsorption of permeants. The DST describes sorption either by simple dissolution leading to mobile and freely diffusible molecules or by an adsorption process producing non-mobile molecules that do not participate in the diffusion process. DST is a quantitative description of penetrant solution and diffusion in microheterogeneous media. This theory postulates that two concurrent modes of sorption are operative in a microheterogeneous medium. The biological tissues, e.g. skin is a heterogeneous structure comprises of keratinaceous cells and lipids, etc. In the cases of the drugs that have low skin permeability, adsorption of drugs in the skin slows down the steady-state permeation. Hence, the total concentration of the drug in the skin can be expressed aswhere,

where,

C = Total drug concentration

C = Mobile drug concentration

= Immobile drug concentration

k = Partition coefficient

C = External concentration

b = Langmuir infinity

C = Saturation constants

Non-linear sorption isotherms can be decomposed into a linear part that accounts for normal dissolution and a non-linear Langmuir-type part that accounts for immobilization of penetrant molecules at fixed sites within the medium [77]. By applying Langmuir equation [66], the concept of diffusivity and flux, the ultimate equation achieved waswhere,

l = Skin thickness

D = Drug diffusivity in absence of drug immobilization.

Chandrasekaran and co-workers [78] explained the skin steady-state flux of an aqueous solution of scopolamine. In vitro flux depicted the linear relationship of concentration with time with a flux of 26 µg/cm2/h at a concentration of 64 µg/ml. The steady-state diffusivity was said to be independent of concentration of the drug, however, lag time was found to increase with an increase in drug concentration. The results of the study were found to be in consonance with the theory.

The investigation carried out by Kubota and his coresearchers also revealed that Timolol permeation could be explained by the dual sorption model [67]. However, it could not be proved that the behavior of finite dose percutaneous permeation is inherent to dual sorption model or not.

A dynamic mathematical model was developed on the basis of bilayer skin i.e. two compartment body model to predict drug plasma concentration after transdermal delivery [79]. The model revealed that effect of metabolism reaction in the viable skin, drug binding, a reservoir in SC, solubility and diffusivity of drug in the skin could be simulated.

Another study revealed that after vehicle removal, dual sorption model still persists. In this case, drug flux versus time depicted straight line graph on a semilog scale as in the linear model, however, half-life was prolonged thereafter. The initial half-life after vehicle removal was found to be longer. This observation confirmed that dual sorption model explains non-linear kinetic behavior of drug [69] (Table 4).

Models based on metabolic activity in the living epidermis and dermis

Several models have been developed based on metabolic activity of drugs in the viable epidermis and dermis. Investigators utilized a computational approach for modeling simultaneously diffusion and metabolism in biological membranes [80]. This approach has a virtually unlimited flexibility in the variety and complexity of metabolic schemes that can be modeled. Furthermore, non-uniform enzyme distributions, as well as composite membranes with any number of layers, can also be handled.

It is expected that this physical modeling approach in combination with other techniques might be a useful tool in the investigation and optimization of drug delivery systems.

Two-dimensional models

In some cases, drug transport through a membrane is considered to be in more than one-dimension. For a two-dimensional model, membrane consists of a biphasic SC and an epidermal/dermal compartment. Several two-dimensional diffusion models have been proposed by different researchers involving different levels of complexity (Table 5). Hence the calculation/prediction of parameters became a more difficult task. Some researchers proposed one of the simplest models considering homogenous SC [81] whereas, other considered the SC as bipasic [85], [86].

Table 5

Applications of two-dimensional macroscopic model.

Reference	Experimental work	Inference
Manitz et al. [81]	Developed a two-dimensional multilayered model representing a cross section of human skinConsidered several penetrating substances formulated within a vehicle for modeling the case of an applied drug and some penetration modifiers (enhancers and reducers, respectively)A coupling via concentration-dependent diffusivities between the diffusion equations of the involved substances was used to model the dependencies between them.	Model equations were solved by exploiting a suitable numerical discretization method.The concentration profile and the position of the penetration boundary were predicted depending on time.
George et al. [82], [83]	Developed a two-dimensional percutaneous absorption modelFeatured a two-dimensional model for non-linear fast adsorption process that is based on the dual-sorption mechanismIncluded a linear second-order parabolic equation with appropriate initial conditions and boundary conditions.Parabolic equations were solved by a finite-difference method which maintains second-order accuracy in space along the boundary.Solution of the parabolic equation gave the concentration of the drug in the skin at a given time	Improved the much simpler one-dimensional models, by giving information of the time-course of drug kinetics in the skin, where the ointment is not directly applied, assuming that the process can be described by a simple homogeneous membrane
Rim et al. [84]	Used a finite element method to simulate two-dimensional drug diffusion from a finite drug reservoir into the skin.Numerical formulation was based on a reaction-diffusion system describing multicomponent nonlinear diffusion which took into account the coupling effects between the different components.	Implemented finite element framework is suitable for modeling of both linear diffusion and nonlinear diffusion with single and multicomponent transport in heterogeneous media where the diffusivities and partition coefficients could vary in each subregion.
Hansen et al. [85]	Developed an in-silico two-dimensional diffusion model consisting of a biphasic SC and a homogenous epidermal/dermal compartmentAn important alternative to in-vivo and in-vitro investigations was considered.Used Fick's second law of diffusion for describing skin transport	Methods and data to measure the partition and diffusion coefficients at various steps during skin transport on the basis of skin anatomical heterogeneity were determined.
Naegel et al. [86]	Developed an in-silico two-dimensional diffusion model consisting of a biphasic SC and a homogenous epidermal/dermal compartmentStudied about the critical parameters during skin transport by comparing experimental and simulated data	Showed how the steady-state flux through SC and interaction of diffusion and partitioning between lipid and corneocytes are related with each other

Microscopic diffusion models

Diffusion models governing transport through a membrane at sub-cellular and cellular level are known as microscopic models. These models can be of two-dimensional or three dimensional geometries. In case of microscopic diffusion models, lipophilic and hydrophilic regimes of SC are represented distinctly. Fortunately, because of improved microscopic methods, it is now possible to collect morphological data at the subcellular level. Most of the macroscopic models rely on certain effective diffusion parameters which are further derived from microscopic models. The incorporation of brick-and-mortar models in structure of SC helped in distinguishing corneocytes (bricks) and the lipids (mortar) that are usually arranged in an idealized geometry. Two-dimensional models involve low computational expense and have the potential of being broadly applicable across diverse chemical classes i.e. hydrophilic and lipophilic substances. However, three-dimensional geometries are closer to the physical reality. They potentially describe the geometrical configurations of the corneocytes and the lipid bilayers as well as the transport processes in the SC in detail, e.g. in combination with swelling. However, these models are computationally expensive. Applications of two-dimensional and three-dimensional microscopic diffusion models are summarized in Table 6, Table 7.

Table 6

Applications of two-dimensional microscopic.

References	Experimental work	Inference
Two-dimensional Microscopic Models
Charalambop-oulou et al. [87]	The most appropriate geometry of porcine SC's lipid and protein phases in a “brick-and-mortar” configuration were determined and further quantified with the barrier properties (diffusivity) of the SC model structures.Spectra of porcine SC are mathematically transformed in order to obtain the corresponding auto-correlation function (ACF).SC structures which matched to this experimentally determined ACF were then produced based on the “brick-and-mortar” configuration.An appropriate numerical method was used to find the effective diffusivity through model domains.	Quantitative determination of the most appropriate geometry of porcine SC's lipid and protein phases in a “brick-and-mortar” configuration and correlate with the diffusivity of the SC model structures.Proposed model could be used for the reconstruction of the intact structure of SC even after hydration.Model was flexible in terms of geometric and other physicochemical parameters.
Frederick [88]	Presented a new mathematical model for permeability of chemicals in aqueous vehicle through skinDiffusion was represented by its fundamental molecular mechanism, i.e. random thermal motion.Diffusion was modeled as a two-dimensional random walk through the biphasic (lipid and corneocyte) SC.Effective diffusivity and effective path length were two main concepts for the prediction of steady-state skin permeability coefficients.	Showed excellent results by regressing four-parameter algebraic model against the Flynn data base.Diffusion phenomena in a morphologically realistic SC structure could be calculated.Provided insight into the fact that SC diffusivity and effective path lengths contribute to overall skin permeabilityHighly complexed computational part and unavailability of limitations to the errors were two major drawbacks.
Frasch and Barbero [89]	Presented finite element (FE) solutions of diffusion through two-dimensional representations of the SC lipid pathwayAssumed that diffusion occurs only within the SC lipids (isotropic nature)	The steady state flux and lag time are solved and compared with the corresponding values for a homogeneous membrane having same thickness consisting of same amount of lipid.Diffusion properties of the SC lipid pathway were correlated with SC geometry.
Barbero and Frasch [90]	Used a finite element method to model diffusion Partitioning was accommodated via change of variables technique.Technique was validated by comparison of model results with analytical solutions of steady-state flux, transient concentration profiles, and time lag for diffusion in laminates.	Diffusion is solved in a two-dimensional “brick and mortar” geometry using partition coefficient between corneocyte and lipid.Results were compared with the diffusion in multiple laminates to examine effects of the partition coefficient.
Barbero and Frasch [91]	Compiled a database of permeability and lag time measurements of hydrophilic compounds from the literature.Used brick-and-mortar geometry to model transcellular and lateral lipid diffusion pathways.	Measured lag times were well correlated (P < 0.0001) with the compound's octanol–water partition coefficient, supporting the role of aqueous–lipid partition mechanism in the permeation of hydrophilic compounds.
Kushner et al. [92]	Developed a new theoretical model – the Two-Tortuosity Model to account for the effect of branched, parallel transport pathways in the intercellular domain of the SC on the passive transdermal transport of hydrophobic permeants.FEMLAB (a finite element software package) was used to validate the model with simulated SC diffusion experiments.	Calculated the vehicle-bilayer partition coefficient and the lipid bilayer diffusion coefficient in untreated human SC.This method was comparable to, and simpler than, previous methods, in which SC permeation experiments were combined with octanol–water partition experiments, or with SC solute release experiments, to evaluate the vehicle–bilayer partition coefficient and the lipid bilayer diffusion coefficient
Wang et al. [93]	Presented a two-dimensional microscopic transport model of the SC named as Wang–Kasting–Nitsche model.Incorporated corneocytes of varying hydration and permeability embedded in an anisotropic lipid matrix.	Dimensionless permeability was calculated as a function of two dimensionless parameters.
Wang et al. [94]	Explicitly considered the sub-structure of the lipid bilayer. Lipid bilayer was treated as a homogeneous medium (with an isotropic diffusion coefficient).	Showed how the model can be related to physicochemical parameters.
Hansen et al. and Naegel et al. [85], [86]	Developed a mathematical model of drug permeation through the SC (brick-and-mortar) and viable epidermis/dermis.Homogeneous diffusivity was maintained in all phases.	Determined the partition coefficients and diffusion coefficients experimentally or derived consistently with the model.Apparent corneocyte diffusivity was estimated by using the apparent SC – and lipid–diffusion coefficients as well as corneocyte–lipid partition coefficients.Quality of the model was evaluated by comparing concentration–SC depth-profiles of the experiment with those of the simulation.

Table 7

Applications of three-dimensional microscopic.

Reference	Experimental work	Inference
Models based on Cuboidal Geometry
Wagner [95]	Cuboid Modeller software tool had been used for the purpose of three-dimensional digital reconstruction of the SC by setting up the model area consisting of hexahedrons with an appearance based on the brick-wall model.Three-dimensional digital SC model areas comprising two homogeneous phases, the lipids and the corneocytes, with variable geometrical parameters were created by Cuboid Modeller.Size of the model area, the size and thickness of the corneocytes, the overlapping degree of the corneocytes and the lipid share were chosen as independent descriptors.	A clearly distinguishable influence of the SC magnitudes on the effective rate of diffusion was established by simulations with different magnitude settings of the SC.Inferred that effective diffusion coefficient was almost proportional to the thickness ratio between the corneocyte and the lipid channel and decreased as the overlapping degree of the corneocytes increasedDistribution coefficient had no effect on the effective diffusion coefficient and the steady state although it considerably influenced the lag time and the diffusion pattern in the non-steady state.The model was observed to improve the basic understanding of the diffusion processes in the SC during the non steady state and the steady state.
Goodeyer et al. [96]	Simulated the steady-state permeation through membranes filled with impermeable flakes numerically using a three-dimensional finite element solver.Compared these numerical results to other published numerical results and the equations that are commonly used to estimate the increased membrane resistance	Inter-layer shoot-through routes can significantly reduced the barrier effect of square-shaped flakes.Introduced the idea of the minimized domain method in which simulations of a single flake layer was combined with a double-half layer calculation to estimate the effective diffusion of any number of layers which inferred that the process can be reversed and the layer-independent effective diffusion can be calculated knowing the effective diffusion for a number of layers of flakes
Marquez-Lago et al. [97]	Simulated water transport and drug diffusion using a three-dimensional porous media modelThe variables that enhanced the complexity of the model were excluded, such as temperature, or variables that arise from the introduction of chemical drugs that would alter its physical structure.	Numerical simulations showed that diffusion takes place through the SC regardless of the direction and magnitude of the fluid pressure gradient.Supported the possibility for designing arbitrary drugs capable of diffusing through the skin whose time-delivery is solely restricted by their diffusion and solubility properties
Models based on Tetrakaidekahedral geometries
Feuchter et al. [98]	Developed a three-dimensional geometry model with TKD for the biphasic model SCTKD were stacked one upon the other and side by side without gaps in a dense packing and with minimal area for all required interfaces.A single corneocyte cell is embedded into a lipid channel yields a TKD base cell T and prepared a cell ensemble which was the agglomerate of several of these cells and these ensembles were horizontally stacked in layers forms a complete diffusional membrane.	Results numerically solved the non-steady-state problem of drug diffusion within a three-dimensional, biphasic model SC-membrane, having homogeneous lipid and corneocyte phases with a multigrid method.
Naegel et al. [99]	Studied the influence of cell geometry on the permeability of the membraneDeveloped a diffusion model based on using numerical simulationCompared the three different geometric concepts (ribbon, cuboid and tetrakaidekahedral type) in two and three space dimensionsGeometric parameters were estimated based on the literature data on human SC morphology.Permeabilities were computed numerically.	The shape of the corneocyte affects the membrane permeability.Results showed that standard brick-and-mortar type models (ribbon, two-dimensional) yielded the lowest permeability values, but exhibited an oversimplified and unrealistic geometry and the model based on TKD (three-dimensional) showed permeabilities that were almost twice larger, but the cell shape is closer to what was observed in reality.
Muha et al. [100]	Evaluated the effective diffusivity in SC membranesTransdermal diffusion was modeled using homogenization theory.Studied the influence of overlapping and swelling of the cells	The method of asymptotic expansion could be used to homogenize membranes consisting of TKD-shaped cells and in calculating the effective diffusivity.Numerical results confirmed that the transversal and lateral diffusivity can be described uniformly with different coefficients.Showed that the efficiency of the barrier increased with increasing overlap and decreased with swelling of the corneocyte cells

Heisig and co-workers [101] presented one of the first two-dimensional diffusion models of SC based on numerical methods. They used a finite volume method for a brick-and-mortar SC geometry which enabled a numerical solution of time-dependent drug concentration in the corneocyte as well as in the lipid phases. Hence, the lag time and drug permeability could be calculated. This model demonstrated that the barrier property of the model membrane depends on three factors; relative phase permeability, corneocyte alignment, and corneocyte-lipid partition coefficient. A representation of a brick-and-mortar geometry is shown in Fig. 3.

Fig. 3

Representation of a brick-and-mortar model (two-dimensional) of SC.

The geometry is defined by the corneocyte width (w), the corneocyte height (h), the lipid channel thickness (d), number of layers (N) and the horizontal cell overlap (ω)where s denotes the horizontal shift.

ω is used to identify the relative factor by which every other cell layer is shifted horizontally. The model helped to predict the diffusion pathways through SC and to identify the diffusivity parameters.

The fit of this model was better than that of the Potts-Guy equation; however, estimation of lag time and skin permeation profiles improved. Wang et al. discussed and confirmed that SC transport pathway is dominated by intercellular hydration and permeability regardless of lipophillicity of components embedding in an isotropic lipid matrix [93], [94]. Nitsche and Kasting developed a 2D microscopic model for viable epidermis layer of human skin on the basis of kerationocytes including tight junctions with adjustable permeabilites to harmonize the transport of solutes (water, l-glucose and hydrocortisone) with lower permeability. The results were represented in terms of effective diffusion coefficient and partition coefficient. However, this theory of Nitsche and Kasting had the limitation as it only defined the basic passive transport properties of viable epidermis of permeation of small solutes. The applicability of such models on large molecules or molecules with high protein binding may be controversial [102].

Three-dimensional models

Models based on cuboid geometry

Unlike previous studies that were limited to consider the two-dimensional cross-section of the SC, homogenization allows the calculation of the three-dimensional effective diffusivity tensor. The modification of the brick-and-mortar models from two-dimensional to three-dimensional was proposed [103], [104]. Rim and co-researchers used the method of asymptotic expansions to formulate the homogenized equation, the effective mass coefficient and effective diffusivity tensor of the three-dimensional SC structure based on a representative unit cell. They performed it for impermeable corneocytes [103] and for permeable corneocytes [104]. Rim and co-workers [103] assumed that the corneocytes were impermeable and studied the effects of different geometric arrangements of the corneocytes within the lipid matrix on the continuum diffusion coefficients. The results indicated that the thin, elongated shape of the corneocytes was ideally suited in lowering the effective diffusivity of the SC. However, Rim et al. considered corneocytes as permeable [104]. The results indicated that the validity of the impermeable corneocyte assumption is largely affected by the solute lipophilicity. Hence, it is demonstrated that hydrophilic solutes with low diffusivities in the lipid matrix exhibit significant mass transport values when the corneocytes are permeable.

Models based on hexagonal prism geometry

Although the cross sections in micrographs of the cuboid models are similar to experimentally observed cross sections, the top view of micrographs depicts a hexagonal cell pattern. A second class modeling of corneocytes as impermeable hexagonal cells was proposed [105]. A Problem Solving Environment (PSE) was developed with connections to remote distributed Grid processes. The Grid is an infrastructure of computing and communications technologies built on the web and on tools for collaborative on-line working [106]. The Grid simulation is a parallel process that allows steering of individual or multiple runs of the core computation of chemical diffusion through the SC. This work demonstrated the feasibility of an intensive finite element solution code run through a PSE. The major drawback of this is the need of efficient handling for visualization of large numbers of data streams.

Models based on tetrakaidekahedral geometries

SC is believed to play the significant role in barrier function of skin. Geometry of cells also influences the membrane permeability. These different geometric concepts were taken into consideration initially, that comprised of ribbon, cuboid and tetrakaidekahedra (TKD) type in two or three dimensions. However, the corneocyte cells have been described by tetrakaidekahedra (TKD) shape. Its geometry uniquely depends upon its height, edge length and the width. Kelvin has described that the clusters of TKD provide a space-filling arrangement with minimum internal surface [107].

Steady-state models

Fick's first law of diffusion described in vitro percutaneous permeation experiments where the steady state is eventually attained after exposing one side of the skin with the solute [108]. It simply relates the amount of solute crossing the skin membrane in a given time. The equation is given as follows:where,

Q = Amount of solute

A = Area of skin membrane

T = Time period for which skin was exposed to the solute

= Concentration gradient constant across the two interior surfaces of the skin

D = Diffusion coefficient

h = Path length

In Equation (27), it is assumed that SC behaves like a pseudo-homogenous membrane which infers that barrier properties do not vary with time. When the solute concentration is maintained at a constant value, the drug flux per unit area from the solute side to the other side, through the skin, approaches the steady-state flux, Jss which is given by:

Now, the maximum flux, , will be acquired when maximum solubility S of the solute in SC will be achieved. Equation (28) will then become

Generally, it is more convenient to express concentrations in terms of the solute concentration in the vehicle (C) and a partition coefficient, K, of the solute between the skin and the vehicle as follows:

These principles are used to demonstrate the predictive value of mathematical models like of skin permeability in defining the absorption of therapeutic and toxic compounds through the skin. Steady-state models are of two types which are described below and their applications are described in Table 8.

Table 8

Applications of steady-state models.

Reference	Experimental work	Inference
Quantitative- Structure Permeability Relationships (QSPRs)
Liou et al. [109]	Developed an empirical model of the QSPR for the transdermal delivery of non-steroidal anti-inflammatory drugs (NSAIDs) to predict the permeability coefficients.Thirteen model NSAIDs were usedIn vitro permeation studies were carried out through the skin of hairless mice. The biological parameters i.e. transepidermal water loss, hydration content, lipid content, resonance running time, and elasticity were measured.	Showed that the solubility parameter might be a more-appropriate drug parameter for predicting the skin permeability of NSAIDs for transdermal deliveryAn empirical model of QSPR to predict permeability coefficient based on the hydrophilicity of the model drugs could be statistically improved with two criteria i.e enhanced solubility and biological parameters of the skin.
Neely et al. [110]	Developed a reliable model integrated nonlinear, QSPR models, genetic algorithms, and neural networksA permeation coefficient data set consisting of 333 data points for 258 molecules was established from literature and the data were added to the extensive thermophysical database.From these data, permeation values for 160 molecular structures were seemed suitable for the modeling.	Established descriptors were employed and new descriptors were built for the aid of development of a reliable QSPR model for the permeation coefficient.An accurate nonlinear QSPR model has been developed for the permeation coefficient and the variability present in the developed model was comparable to the level of variability expected in the experimental data.
Cronin et al. [111]	Permeability coefficient values of 114 compounds across excised human skin in vitro were noted from literature.Forty-seven descriptors were calculated by using relevant physicochemical parameters of the compounds.QSPRs were developed using least-squares regression analysis.Compared with a QSPR based on penetration across a synthetic polydimethylsiloxane membrane	Results indicated that the percutaneous absorption is mediated by the hydrophobicity and the molecular size of the penetrant.Comparison of QSPRs based on penetration across a synthetic (polydimethylsiloxane) membrane inferred that the mechanisms of drug flux across polydimethylsiloxane membranes and excised human skin are significantly different.
Patel et al. [112]	Developed QSARs for the skin permeability coefficients of 158 compounds through excised human skin in vitroMany compounds, including hydrocortisone derivatives, were removed from the dataset as the permeability data for these compounds was showed as an error.Hydrophobicity, molecular size, and hydrogen bonding were selected as descriptors	Descriptors provided an excellent fit to the data (r2 = 0.90). Observed that a good QSARs could be developed utilizing hydrophobicity and molecular size only
Chang et al. [113]	Developed a QSAR for permeability coefficient prediction and for characterization of the physicochemical properties involved in transdermal transport of chemical at a molecular levelFor developing this QSAR, 158 chemical substances with known permeability coefficient were selected.The electrostatic interactions between electric quadrupoles of van der Waals forces, octanol–water partitioning of solute, similarity in antineoplastic property and frequency of carbon–nitrogen bonding at a constant topological distance were selected as four molecular descriptors.	The regression was found to be 0.828 with the four-descriptor multiple linear regression model fit and a mean percentage error of 18.8% was observed which suggested that this QSAR can be used as an alternative source of permeability coefficient information and a potential tool for dermal hazard characterization when experimentally determined partition coefficients are not readily available.
Kang et al. [114]	The human skin penetration effect of 49 terpenes and terpenoids on in vitro permeability coefficients of haloperidol was evaluated.The permeability coefficient of the drug to the lipophilicity, molecular weight, boiling point, the terpene type and the functional group of each enhancer were linked to each other using a first-order multiple linear regression model.The QSAR model was derived from the data generated by using standardized experimental protocols.	Study suggested that terpenes which possess one or combinations of the properties like level of hydrophobicity, phase (liquid state), appearance of specific functional groups (ester or aldehyde but not acid), and chemical types (not a triterpene or tetraterpene) might be better enhancers for drug permeation through skin.Useful in designing of new terpene enhancers and for preliminary screening of terpenes as penetration enhancersQSAR models allowed the prediction of human skin penetration effect for other terpenes and terpenoid compounds for drugs having physicochemical properties similar to heloperidol without the need to conduct in vitro experiments with human skin samples which are not easily available.
Lai and Roberts [115]	Examined the extent of the ionic mobility–pore model which could be used to describe epidermal iontophoretic structure–permeability relationshipsIncluded iontophoretic transport, solute size, solute mobility, total current applied, presence of extraneous ions (determined by conductivities of both donor and receptor solutions), permselectivity of the epidermis, as well as a solute pore interaction as determinants	All of the determinants together provided an excellent regression for iontophoretic permeability.The pore radius estimated from analysis of cations, anions and uncharged cation of the ionic mobility–pore model to the transport of solutes ranged from 6.8 to 17 A°.
Polat et al. [116]	Investigated the effect of sodium lauryl sulfate (SLS) on skin structural deviations along with low-frequency sonophoresis (LFS)Pig full-thickness skin and pig split-thickness skin could treated with LFS/SLS and LFS and then analyzed using aqueous porous pathway model to quantify skin deviations by observing changes in skin pore radius and porosity-to-tortuosity ratio.Skin electrical resistivity was also analyzed to find out the effect of SLS on reproducibility and predictability of skin perturbation.	Showed that LFS/SLS-treated pig-full thickness skin, LFS/SLS-treated split-thickness skin and LFS-treated full-thickness skin exhibited similar skin changes but LFS-treated split-thickness skin exhibited significantly higher skin changes.Evaluation of porosity-to-tortuosity ratio values suggested that both types of the skin when treated with LFS/SLS have similar transport pathways but when they are treated with LFS alone, they possess a low porosity-to-tortuosity values.It was concluded that the simultaneous use of SLS and LFS not only results in synergistic enhancement, but also in more predictable and reproducible skin changes.

Quantitative structure-permeation relationship (QSPR) models

A QSPR model, also known as quantitative structure–activity relationship (QSAR) describes a mathematical relationship between structural attributes and a property of the set of chemicals [117]. The use of such mathematical relationships to predict the target property of interest for a variety of chemicals prior to or in lieu of expensive and labor intensive experimental measurements has naturally been very enticing. The potential promise of using QSPR models for screening of chemical databases or virtual libraries before their synthesis appears equally attractive to chemical manufacturers, pharmaceutical companies and government agencies, particularly in times of shrinking resources. Many different methods embedded in a variety of computer software packages are available for this purpose; the choice of software and model is generally guided by the answers or solutions of the resulting QSPR screen is expected to provide. Some widely used model development methods include multiple linear regression (MLR), partial least squares (PLS), artificial neural networks (ANNs), and k-Nearest Neighbour (kNN) methods. For reliable predictions, QSPR models should be statistically significant and robust and should be validated and application limitations should be defined.

The main focus of QSPRs has been the assessment of a permeability coefficient, which is defined as the steady-state flux of chemical across the skin normalized by the concentration gradient [14],where,

k = Permeability coefficient.

Flynn published a large data of skin permeability coefficients of a number of drugs from aqueous solutions [118]. Potts and Guy [6] then took advantage of this data to generate a QSPR that is now the most cited and applied QSPR model for predicting skin permeability:where,

k is in cm/s

MW = molecular weight of penetrant/drug ranging from 18–750

log p values are from −3 to +6.

The incorporation of some modifications in terms of the compound's molecular volume, polarizability, hydrogen bond donor and acceptor activities and molar refractivity, Potts and Guy improved the fit of model to more limited datasets [119]. Another QSAR model involving partition coefficient, molecular weight and hydrogen bonding was utilized for correlating the structures and skin permeability of a wide range of compounds through human skin and through hairless mouse skin [120]. The correlations obtained were found to be dependent on biological systems as well as on the physicochemical properties of vehicle used. Few researchers utilized parameters of free energy of solute transfer to lipid related to molecular structure instead of lipophillicity [119], [121].

Although most of the QSPR models are used to predict permeability coefficient (k) but is practically more a relevant parameter. It can be estimated using experimental values of the aqueous saturation concentration along with k calculated from a QSPR model for an aqueous vehicle. Hence, based on these assumptions QSPR model for , rather than k, was also developed [122]. Interestingly, here log P seemed to be a less significant parameter and MW alone was sufficient to describe the bulk of the chemical specific variation in the data. A well known “500 dalton” rule states that the MW of a compound must be under 500 daltons so as to get permeate through the skin. Thus, compounds having MW greater than 500 daltons cannot pass through the corneal layer and compounds with smaller MW can easily surpass transcutaneously.

The development of rational and improved QSPR models is continuously attracting the attention of researchers. The strategies employed in QSPR are quite novel, and QSPR may become a new and better tool to solve the problem. Such approaches include chemical structure-based approaches, for example, fragment descriptor and neural network-based modeling approaches using large databases of theoretical molecular descriptors [123], [124], [125], ensemble modeling using nearest neighbor theories [126], topostructural, topochemical, shape and/or quantum mechanical indices [127] and Gaussian process models [128]. Such approaches, though novel, have yet to find significant, “real-world” application.

Although QSPRs have many benefits, they have certain limitations as well. Most importantly, they cannot be used when formulation components modulate the barrier properties of the skin. Extrapolation of the predictions from simple aqueous solutions to complex multicomponent and/or multiphasic formulations is ambiguous and may be unachievable in some cases. However, QSPR modeling along with computer aided molecular design technique may demonstrate designing or generation of potential chemical penetration enhancers beyond the chemical structure available commercially. Further, the established QSPR models may allow the use of structural modifications in existing permeants/enhancers in predicting existing permeation or structure modifications in drugs that may lead to increased permeation without the need to conduct in vitro experiments.

Porous pathway models

Generally, most of the equations which are based on permeation through lipids adequately describe the permeation of lipophilic drugs across SC but their applicability to hydrophilic drugs is not accurate. Many attempts have been made to develop models to describe the transport of hydrophilic molecules. Appendages (hair follicles and sweat ducts) are a good pathway for permeation of hydrophilic solutes. The average density of hair follicles in human skin is about 50–100 cm−2 and the area fraction occupied by the follicles is about ~10−3 [129], [130]. The major part of the hair follicle area is occupied by the hair shaft. The sweat glands occupy an area fraction of about 10−4. If the diffusion coefficient of a solute through appendages could be ~10−6 cm2/s with an available area fraction of appendages of 10−4 cm2, then the contribution of appendages to skin permeability of the order of 10−6 cm/h [131]. However, many hydrophilic solutes permeate skin at a much faster rate [132]. Tezel and co-workers elucidated the incorporation of lipophilic pathway into the process pathway model; in absence and in presence of ultrasound in order to explain the permeation through skin [133].

A model termed as a porous pathway model was introduced to find the permeability of hydrophilic solutes [134]. The imperfection in the SC lipid bilayers is the rationale behind the existence of “pores”. These imperfections may manifest themselves as separation of lattice vacancies, multimolecular voids due to missing lipids, grain boundaries, or defects created by steric constraints placed by the corneocytes on intercellular lipid bilayers. It might be possible that these imperfections offer low energy channels for diffusion of solutes across poorly permeable intercellular lipid bilayers. The general expression for permeability based on the porous pathway, kp of a hydrophilic permeant is given by:where,

ε = Porosity

τ = Tortuosity

= Thickness of the membrane

= Diffusion coefficient of the permeant in the liquid-filled pores within the membrane.

According to the hindered transport theory, is a function of both the permeant and the membrane characteristics. expressed as a product of the permeant diffusion coefficient at infinite dilution and permeant diffusion hindrance factor which depends on the pore size [135]. The generation of porous pathway through low frequency ultrasound perturbed the structure of skin in a reproducible and predicted manner more significantly than the use of SLS [136]. All these investigations supported the fact that creation of pores in the structure of skin could enhance the rate of permeation. The pores can be generated as a result of use of chemical enhancers or physical techniques.

Transient models

Previously discussed models are based on steady state permeability across a membrane beyond the diffusion lag time. There are few models that have described the relationship of time and skin penetration.

Basic models

Transient drug diffusion across the SC is given by Fick's second law as given below:where,

C = Concentration of the permeating solute at time t

x = Depth of the skin to which permeant has reached.

It is assumed that the SC behaves like a pseudo-homogenous membrane and the diffusion coefficient and partition coefficient do not vary with time. To solve this equation, the starting concentration within the SC as well as the conditions for concentration or flux at the boundaries of the SC (i.e., at the outermost and inner most surfaces of the SC) must be specified.

Compartment models

Mathematical models are a collection of mathematical quantities, operations and relations together with their definitions and they must be realistic and practical. The compartment is the traditional and most widely used approach to characterize the drug pharmacokinetically. Hence, compartmental models are also known as pharmacokinetic (PK) models. These models are often used to study the fate of chemicals entering and leaving the body. These models simply interpolate the experimental data and use empirical formulas to estimate the drug concentration with time. These models help in characterizing the behavior of the drug in the patient, explain the drug interactions, determine the influence of altered physiology or disease state on drug ADME and help in calculating the optimum dosage for patients. Compartmental model, in combination with pharmacodynamic models, are potentially useful tools for risk assessments and predictions of transdermal drug delivery. These models treat the skin and also the body as one or several well-stirred compartments of uniform concentration that act as reservoirs of chemical storage and transfer occurs within the compartments depicted by the first order rate constant expressions.

A differential mass balance equation of chemicals in the one-compartment skin-layer gives the equation as follows:where,

= Position-averaged drug concentration in the skin layer

= Volume of the skin layer

k = Rate constants describing drug transfer between the vehicle, skin and blood compartments (j = 1, −1, 2 and −2).

The two-compartment model of skin is given by:where,

= Position-averaged drug concentrations in the SC

= Position-averaged drug concentrations in the viable epidermis

The benefit of compartment models is that the mathematical solution of even complex exposure situations combined with variable distribution and metabolism in and elimination from the body as well as variations in blood flow to the skin is represented by first-order differential equations that are easily solved by a number of standard software packages. These models are also simple enough to conduct probabilistic calculations, allowing an assessment of the effects of variation in the system parameter. Many researchers have used compartment models to evaluate the data (Table 9).

Table 9

Applications of compartmental model.

References	Experimental work	Inference
Wallace and Barnett [137]	A compartmental model was developed to describe the penetration of a drug from a topically applied vehicle through the skin.Data for in vitro penetration of methotrexate through hairless mouse skin from vehicles having varying pH from 3.5 to 6.5326 was computed and fitted to estimate model parameters	Observed the existence of parallel penetration pathways exist by comparing lag time and the exponential coefficientDrug penetration increased through the shunt pathway as the pH and ionization of vehicle increased.Flux through the shunt pathway was maximum at pH 6.5.
Kubota and Maibach [138]	A compartmental model was proposed to generate finite-dose percutaneous permeation pharmacokinetics.Compared the lag-time and steady-state flux predicted by these models with those predicted by the diffusion model	Found that three statistical parameters of the compartmental model viz. the mean residence time of the drug in vehicle and that in skin and the variance of the residence time in the vehicle are same to that of diffusion modelSame lag-times and steady-state fluxes were given by both of the models.
Nugroho et al. [139]	Developed a number of compartmental models to describe the transdermal iontophoretic transport of drugs in vitro quantitatively.Two different compartment models were described for in vitro transport during iontophoresis and one compartment model for in vitro transport in post-iontophoretic period.Skin was used as intermediate compounds and all models were based on mass transfer from donor to receptor compartment.	After computer-fitting of all parameters i.e. lag-time, steady state flux during iontophoresis, skin release rate constant, the first-order rate constant of the iontophoretic driving force from the skin to the receptor compartment and passive flux in the post-iontophoretic period found during iontophoretic and post-iontophoretic periodEffect of iontophoretic driving force on mass transfer from donor to receptor compartment was not found to be significant.
Nugroho et al. [140]	Models were based on in vitro compartmental models described by Nugroho et al.Developed two different models based on the mode of input of drug into systemic circulation (One model was based on constant input and other was on time-variant input).Four parameters were taken i.e. lag time, steady-state flux during iontophoresis, skin release rate constant and passive flux in post-iontophoretic period.Elimination was described by a one-compartment and two-compartment model.	Observed that the time-variant model was more appropriate than the constant input models.Provided the basis to find in vivo-in vitro correlation of drug transport by iontophoresis
Saurabh and Chakraborty [141]	Developed a three-compartmental mathematical model of an implantable polymer membrane-based drug release device in which they used different voltage across the electro-active membrane for delivery and uptake of anionic anti-tumor drugsConsidered the drug-reservoir, the polymer membrane and the diseased tissue as three compartments and used Poisson–Boltzmann, Nernst–Planck and Diffusion–reaction models for these compartments, respectivelyUsed Laplace transforms and residue integration to obtain an analytical solution	Quantified the various controlling parameters and analyze the dynamics and efficacy of the delivery method of an anionic anti-tumor drug, Irinotecan-HCL, commercially available as CPT-11The drug concentration in the tumor tissue could be maintained within the therapeutic range by adjusting the applied voltage.

Conclusion

By using mathematical models, permeation data can be obtained without performing experimental work, shunning the ethical issues concerned with the use of animals. Also, tedious and time consuming experimental work can be fended off by using these models. As the number of models is bountiful and each model is very complex, hence, this overview is limited in its ability to give each of the models the credit they may deserve. However, more emphasis is being given on the models that have been used practically by various researchers in recent studies to explain this complex domain.

Mathematical modeling has a momentous possibility to facilitate product development in the future and to help understanding the complex pharmaceutical dosage forms. Due to the new developments in information technology, the accuracy of these models has been steadily increasing and they are becoming easier to employ. Mathematical modeling of drug delivery can be expected to become an integral part of product development like other scientific disciplines. However, major drawback of mathematical modeling is that there will be one general theory which will be applicable to any type of drug delivery system. This drawback can be overcome by developing a broad spectrum of different mathematical models, applicable to specific types of devices differing in geometry, drug and excipient type. Decisive parameters will help in identifying the appropriate model for a specific type of delivery system and type of task e.g., prediction of the effects of formulation parameters or understanding of the underlying drug release mechanisms.