Predicting the transport of chemicals across the skin: Stochastic calculus and machine learning

Highlights

• Stochastic differential equations (SDEs) were applied to model the inherent variability in skin permeability data, addressing limitations of deterministic quantitative structure– permeability relationships (QSPR) models.

• Three stochastic approaches – the Euler–Maruyama algorithm, the Milstein scheme, and the Heston model – were implemented and compared for predicting permeability coefficients.

• Six molecular descriptors (molecular weight, log Kow, Balaban Index, Harary Index, Ramification Index, and Forgotten Index) served as predictive inputs.

• The Euler–Maruyama method demonstrated superior stochastic modeling performance with an R² of 0.81 and MSE of 0.5259.

• Among the machine learning approaches, gradient boost regression achieved the best predictive accuracy (R² = 0.84, MSE = 0.45), outperforming both support vector regression and the Levenberg–Marquardt algorithm.

• Comparative analysis revealed that gradient boost regression and the Euler–Maruyama method offer the most reliable predictions for skin permeability coefficient estimation.

• This study provides a framework for integrating stochastic modeling with machine learning to account for experimental uncertainty in permeability predictions.

INTRODUCTION

Skin permeability coefficient (log Kp) prediction is essential for the design of transdermal medication delivery devices and for risk assessment.^[1] A variety of in silico and in vitro experimental methodologies can be employed to forecast and quantify, respectively, the permeation of chemicals through the skin.^[2] The mathematical tools comprise equations employed to forecast skin permeability based on physicochemical characteristics, including molecular weight (MW) and the octanol-water partition coefficient (logKow), among others.^[2] The prediction of cutaneous absorption is a significant study area in the pharmaceutical and cosmetics industries, facilitating the optimization of transdermal permeation and aiding in hazard and risk assessment.^[3] The primary advantages of theoretical predictions compared to experimental measurements encompass economic efficiency and the mitigation of ethical dilemmas.^[3] Furthermore, the models may facilitate a more profound comprehension of absorption mechanisms.^[3] The health risk that comes from skin contact with harmful chemicals can be measured by finding the substances’ skin penetration coefficient (Kp). Nonetheless, due to inherent limitations of the technique, Kp evaluation is not performed satisfactorily.^[4] Despite the remarkable advancements achieved in recent decades, precise mathematical modeling of skin permeability continues to pose a significant challenge.^[5]

A crucial factor determining a molecule’s capacity to permeate the skin is its permeability coefficient (Kp).^[6] Quantitative structure–permeability relationships (QSPR) generally employ descriptors^[7] such as MW, partition coefficient (log Kow), and topological indices to forecast log Kp.^[8] While QSPR models have enhanced our comprehension of structure–permeability correlations, their deterministic characteristics frequently underestimate the variability observed in actual datasets,^[9,10] which stems from variations in skin morphology, experimental circumstances, and measurement noise.^[11,12] From a therapeutic perspective, it is essential to establish methodologies that can characterize and forecast the uncertainty in transdermal permeation parameters. Uncertainty quantification is receiving heightened attention, and numerous strategies are being devised to address this challenge.^[13] Epidermal permeability, commonly represented by the stratum corneum permeability coefficient, Kp, is a significantly variable experimental parameter.^[5] The measurements are influenced by the skin’s state and the experimental design.^[5] In analogous in vitro tests, Kp-values may fluctuate by several orders of magnitude based on the source of the skin sample, such as whether it is derived from different persons, cadavers, or anatomical regions.^[5] Similarly, experimental variables such as temperature, vehicle, pH, and others also affect Kp-values.^[5]

Classical QSPR models offer deterministic predictions from molecular descriptors, but they often neglect biological variability and measurement uncertainty.^[9,14] The predominant erroneous assumption in quantitative structure-activity relationship (QSAR) and quantitative structure-permeability relationship (as well as QSPR) modeling is that the reported value for any experimental endpoint represents the definitive “true” value for that measurement.^[9] It is assumed that the provided experimental value represents the sample mean, which adequately approximates the population mean (actual value) of all potential measurements.^[9] Consequently, data modeling often fails to facilitate a comprehensive knowledge of population distribution and associated uncertainties.^[9] Uncertainty in physico-chemical data is intrinsic, regardless of whether the data are measured or modeled.^[14]

Differential equations have been utilized by researchers to characterize systems.^[15] Several scholars have employed ordinary or stochastic differential equations (SDEs) to model various dynamical systems.^[16] Stochasticity, or unpredictability, influences mathematical modeling of scientific predictions.^[17] SDEs are mathematical models extensively employed to characterize complex phenomena that are sometimes influenced by noise.^[18] The dynamics of biological systems are predominantly affected by stochastic processes and are susceptible to random external disruptions.^[19] The ramifications of such processes are frequently examined through the formulation and evaluation of stochastic models.^[19] A multitude of scientists have employed SDEs to model diverse events.^[20] SDEs can concurrently encapsulate the established deterministic behaviors of relevant variables while allowing a modeler to account for the unexplained random variables.^[20]

Conventional QSPR models frequently employ deterministic regression based on molecular descriptors, including MW and lipophilicity. The intricate and diverse nature of the skin barrier, along with the random characteristics of molecule diffusion, necessitates the application of SDEs.^[18] The Euler–Maruyama (EM) method is a widely used numerical methodology for modeling SDEs.^[21] We additionally employed the Milstein scheme (MS) for forecasting. The Milstein approach approximates solutions to SDEs by integrating higher-order Itô’s lemma terms, hence enhancing accuracy compared to basic Euler–Maruyama techniques.^[22] MS improves QSPR models of skin permeability by including uncertainty while preserving predicted accuracy.^[13] It offers probabilistic predictions that facilitate decision-making in transdermal medication development and safety evaluations. In the Heston model, volatility is permitted to “diffuse” throughout the entire space governed by a diffusion mechanism.^[23] The Heston model has two correlated SDEs.^[23] The spectrum of Heston’s volatility is continuous. This study project employed the EM algorithm, the MS, and the Heston stochastic technique^[13,24,25] to forecast the permeability coefficient.

In addition to the above-mentioned stochastic models, we also utilized machine learning algorithms (Gradient Boost, Support Vector Regression [SVR], as well as the Levenberg– Marquardt [LM] technique) for prediction. Gradient boosting (GB) methods use numerical optimization for classification and regression by iteratively advancing in the approximate direction of the negative gradient of the loss function.^[26] Gradient boosting strategies typically utilize a “weak learner” (such as a regression tree) as an approximation of the steepest descent path, owing to the difficulties associated with accurately navigating the negative gradient.^[26] The Gradient boosting regression (GBR) is an ensemble model comprising an iterative sequence of increasingly structured tree models, with each successive model learning from the errors of its predecessor.^[27] This machine learning approach employs “boosting” of a collection of weak prediction models, usually decision trees (DTs), to develop a more resilient model. A GBR of M trees is conventionally represented as:

Equation 1

where f_m is the predicted output value, M is the number of DTs, γ_m is a scaling factor adding the contribution of a tree to the model, and h_m is the weak learner.^[27] The forecast for a new observation is obtained by aggregating the predictions of all individual trees that comprise the model.^[28] Gradient boosting utilizes a statistical method to improve traditional DT models.^[29] The objective is to integrate a set of underlying models to develop a unified, robust model. The GB approach progressively constructs new DTs through residual reduction. This iterative method improves estimates by progressively integrating a new tree with a minimized loss function.^[29]

SVR is a widely utilized supervised machine learning method.^[30] It serves as the regression equivalent of the more widely recognized support vector (SV) machine classifiers.^[30] SVR can be employed for non-linear regression (LR) through a technique referred to as the kernel trick.^[31] This regression approach aims to estimate a function f(x) based on the principle of structural risk minimization derived from statistical learning theory.^[32] It improves generalization capability by decreasing the upper limit of risk, hence alleviating overfitting and leading to a diminished generalization error.^[32] The SVR model produces markedly nonlinear functional relationships between the input data and the feature space.^[30] This technique enables comparatively straightforward mathematical formulations and calculations.^[30] The radial basis (RB) function (RBF) kernel function converts the original input space into a feature space with higher dimensions.^[33] The aim is to determine the regression hyperplane created by the SVs

Observations in the calibration set nearest to the hyperplane are referred to as SVs.^[34] The direction of the hyperplane is defined by SVs. Epsilon (ε) quantifies the allowable error in the model’s predictions. An elevated ε permits a wider buffer for accommodating larger prediction errors.^[34] The cost (C) regulates the equilibrium between error reduction and hyperplane margin enhancement, while sigma (σ) can alter the Gaussian basis functions in the RBF kernel. A diminished σ results in a more complex model.^[34] SVR utilizes an error measure characterized by a margin of tolerance (ε-insensitive) to ignore errors that are less than ε.^[35] This makes the fitting less vulnerable. SVR differs from LR using an ε-insensitive loss function instead of the traditional L2-norm loss function, hence improving its outof-sample performance.^[32] The principle of SVR is to identify the best hyperplane that contains the maximum number of data points within the data space.^[33] The effectiveness of an SVR model depends on three essential factors.^[36] The parameter C enables the balance between training error and model complexity. Excessively large values of parameter C will result in overfitting. If C is too small, it may result in underfitting and an increase in training errors.The parameter ε governs the width of the ε-insensitive zone, employed to accommodate the training data. An increased ε-value results in a diminished quantity of selected SVs and yields flatter or less intricate regression estimations. If ε is very big, the separation error is substantial, the number of SVs is low, and vice versa. Nonlinear mapping from the input space to a high-dimensional feature space is determined by the kernel parameter (γ).^[37]

The LM curve-fitting method combines two minimization methods: the Gauss-Newton method and the gradient descent method.^[38] This optimization technique has been used to solve non-linear least squares problems. The identification of a model from a large dataset depends on the parameter optimization method.^[36] Gradient-based optimization methods are commonly utilized for solving nonlinear problems because of their high reliability and computational efficiency.^[36] The LM method is an iterative approach that can be used to compute unknown parameters.^[39] LM includes a damping component that facilitates the dynamic modification of the iterative direction and step size, bridging the Gauss–Newton method with the steepest gradient technique.^[36] By carefully selecting the damping factor, suitable convergence can be achieved, hence avoiding the constraint of a full rank (rows or columns are all linearly independent) that is required for a Jacobian matrix (a matrix of all first-order partial derivatives of a system).^[36,40]

QSPR models connect some structural properties of molecules (termed molecular descriptors to skin permeability.^[41] A topological index is a number assigned to a molecule based solely on its connectivity pattern or topology without taking into consideration angles or bond lengths.^[42] Topological indices (a subset of molecular descriptors) are important tools in QSPR studies.^[43] A molecular graph represents the structural formula of a chemical compound, where the vertices correspond to atoms and the edges represent atomic bonds.^[44] Molecular graphs are applied in cheminformatics, QSAR, and QSPR.^[45] Numerous topological indices have been proposed and used in various research domains.^[46-49] Topological indices can be categorized into three groups: distance-based indices, degree-based indices, and spectrum-based indices.^[50] A multitude of distinct descriptors (also called invariants) numbering in the hundreds have been utilized in QSAR/QSPR studies.^[51] Degree-based indices encompass the Zagreb index, Randić index, and connective eccentricity index, among others.^[52] Distance-based indices include the Wiener polarity index, Wiener index, Szeged index, Mostar index, Kirchhoff index, ABC index, Harary index, among others.^[53,54] Spectrum-based topological descriptors refer to those computed from eigenvalues of valency-related or distance-related chemical matrices.^[55] The adjacency energy, signless Laplacian energy, the Estrada index, the distance energy, as well as the Laplacian energy are spectrum-based topological indices.^[55]

We employed topological indices as features (independent variables) in the machine learning models for predicting permeability. Topological descriptors are numbers used to analyze various physicochemical properties of molecules.^[56] They are generally derived from chemical graphs. Chemical graphs portray molecules, with atoms depicted as vertices while edges represent chemical bonds.^[57] According to the chemical graph theory, certain molecular properties can be evaluated using the data reflected in their corresponding chemical graphs.^[57] A graph invariant, also known as a topological index, is a function applied to a graph that remains unaffected by the labeling of its vertices.^[57] A molecular graph is a fundamental graph linked to a molecular structure.^[49] Each vertex of a molecular graph represents an atom of the molecule, while the edges denote the atomic bonds.^[49] Certain topological descriptors are utilized by researchers to evaluate drug molecules, facilitating the identification of important pharmacological features.^[49] A molecular graph G=(V, E) consists of the edge set E(G) and the vertex set V(G).^[49] The degree of a vertex v in graph G, denoted as dG(v), is the number of edges incident to it.^[49] The topological index of a graph G is Top(G), and that of a graph H is Top(H). When Top(H) = Top(G), then graph H is considered to be isomorphic to graph G.^[49]

The Balaban index of a graph G is a distance-based topological metric that has been effectively employed in many QSPR and QSAR models.^[58] It has been demonstrated in comparative studies investigating the Balaban index [J(G)] and the Wiener index [W(G)] for alkanes that the J index decreases the degeneracy of W and imparts markedly improved discriminative power.^[58] The Balaban index is an important topological index in chemical and biological research.^[59] It is comparable to the Wiener index, but displays diminished degeneracy and improved discriminative power.^[59] This is especially beneficial for more complicated and larger molecular structures. Degeneracy occurs when the topological metrics do not differentiate between molecules with equal cycle numbers or atom counts. The Balaban index solves this issue more successfully than competing indices, showcasing its relevance. This can be attributed to its high discriminatory ability and low degeneracy. J is especially useful when investigating complex molecules with numerous heteroatoms and bonds. Several researchers have used this index to study numerous properties such as octane ratings, density, boiling point, and refractive index. It has also been employed in the study of supramolecular substituted sulfonamides.Some investigators have used this index to predict the melting and glass transition temperatures of linear macromolecules and to determine the properties of polymers. J has been utilized in biological research that focuses on optimizing physiologically active lead compounds, alongside other chemical descriptors.^[59]

The Harary index [(G)] is a distance-oriented topological index.^[60] (G) was independently proposed in 1993 by Ivanciuc et al.^[61] and Plavšić et al.^[62] In a connected molecular graph, the Harary index is the inverse of the total distance sum between all unordered vertex pairs of G.^[63] The Harary index is a useful measure for evaluating the compactness of a molecule. An elevated Harary index signifies enhanced molecular compactness. The Harary index serves to assess branching in alkanes. In compact structures, the number of short distances increases, leading to a Harary index that surpasses that of extended structures. Similarly, a more complex configuration demonstrates greater short distances, leading to an elevated Harary index for these structures.^[64]

In 2015, Furtula and Gutman introduced a topological index, usually referred to as the Forgotten Index.^[65] The Forgotten topological index is denoted as the summation of the squares of adjacent vertices. This index is regarded as a suitable instrument for examining the medical and biological data of drug molecular structures.^[66] Furtula and Gutman proposed the forgotten index (F-index) to predict the characteristics of octane isomers.^[65] This index was referenced in the same article that presented the Zagreb indices in 1972, albeit it garnered less notice.^[67] After 43 years, Gutman and Furtula noted that the overlooked index is an extraordinarily effective and powerful instrument for ascertaining the characteristics of organic and inorganic compounds.^[67] The F-index is defined as the summation of the cubes of the degrees of the vertices in molecular networks.^[68] The Ramification Index is computed as the summation of all vertex degrees exceeding two, subtracting two from each vertex degree.^[69,70] In this project, we used six different regression approaches. Table 1 shows the comparative analysis of stochastic and machine learning techniques. A working model of this project is shown in Figure 1 and Appendix 1.

Close

Figure 1:

Graphical abstract. MSE: Mean squared error, GBR: Gradient boost regression, E-M: Euler–Maruyama algorithm, SDEs: Stochastic differential equations.

Export to PPT

MODELING FRAMEWORK FOR PERMEABILITY COEFFICIENT PREDICTION

METHODS

RESULTS

Simulation and modeling were conducted with techniques from stochastic calculus. Figure 2 shows a parity plot that compares observed log Kp values with the EM mean predictions. The plot from the EM method shows strong agreement between observed and predicted values across the logKp range (−6 to 0). Most predictions cluster near the identity line, with modest deviation at extreme values (very low permeability compounds). The data suggest the EM approach accurately captures permeability trends. The slight underprediction at extremes can be attributed to the limits of molecular descriptors or insufficient sampling in high-variance regions. An MSE of approximately 0.43 and an R² of about 0.82 confirms the predictive accuracy of the EM model that is comparable to deterministic QSPR models, with the added benefit of variance estimates. By incorporating stochastic variability, the EM method avoids overconfidence in permeability predictions. This is especially critical for risk–benefit analysis in regulatory submissions, where uncertainty quantification strengthens robustness. The EM approach offers several key benefits. Rather than providing only a point estimate, the EM simulation yields a predictive distribution, allowing calculation of credible intervals for log Kp. Second, using a QSPR-based regression as the drift term, the method merges mechanistic insights from molecular descriptors with the flexibility of stochastic calculus. The method can be extended to multi-zone skin models and time-dependent release kinetics. Despite these advantages, some limitations should be recognized. Accurate estimation of theta and sigma is essential; poor calibration can distort both mean predictions and uncertainty estimates.

Close

Figure 2:

Parity plot that compares observed logKp values with the Euler–Maruyama (EM) mean predictions.

Export to PPT

When using the Milstein method, a Random Forest (RF) model was used to provide baseline predictions of logKp. Residual variance defined the diffusion term. Milstein simulations (2000 paths per compound) generated distributions of predicted logKp values, providing both mean predictions and confidence intervals. The Milstein stochastic augmentation preserves nearly all of the baseline model’s predictive accuracy while adding a probabilistic framework that provides confidence intervals and quantifies uncertainty. The very small changes in R² and MSE show that incorporating the Milstein SDE does not degrade the core prediction quality; instead, it enriches the model with meaningful measures of variability that the deterministic baseline cannot supply. Minor deviations at the extremes (e.g., highly negative log Kp) reflect increased difficulty in modeling compounds with very low permeability, where experimental variability is greatest. The plot supports the conclusion that the EM mean retains the accuracy of the baseline regression while offering additional uncertainty quantification. Figure 3 (observed vs. predicted logKp, Milstein mean) compares experimental logKp values on the x-axis to the mean of the Milstein simulations on the y-axis. The dashed 45° line represents perfect agreement. Most points lie close to this line, indicating that the Milstein approach accurately reproduces experimental permeability coefficients. The scatter around the line reflects both biological variability and model error, and is consistent with the 95% confidence intervals derived from the stochastic simulations. Figure 4 (predictive uncertainty) shows a histogram of the widths of the 95% confidence intervals for predicted logKp values. Most intervals are about 1.7–1.9 log units wide, illustrating that the Milstein method provides a stable quantification of predictive uncertainty across different compounds. This information is valuable for risk assessment and for selecting candidates with acceptable variability in skin permeation.

Close

Figure 3:

Observed versus predicted logKp, Milstein mean.

Export to PPT

Close

Figure 4:

Predictive Uncertainty: a histogram of the widths of the 95% confidence intervals (CI) for predicted logKp values.

Export to PPT

The Milstein Sample Paths [Figure 5] show multiple stochastic trajectories of logKp versus time for a representative compound. Each colored line represents one Milstein-simulated path, while the dashed horizontal line marks the experimental logKp. The ensemble of paths captures the stochastic fluctuations around the mean prediction and confirms that the experimental value lies within the simulated envelope. This illustrates how the MS captures the inherent variability of the diffusion process. The Milstein method provides several advantages over simpler stochastic schemes. It accounts for state-dependent noise through the inclusion of the derivative term b’(X), offering improved strong convergence. For skin permeability modeling, this means that variability in logKp predictions scales realistically with molecular properties. The predictive distributions and confidence intervals produced enable more informed decisions in pharmaceutical formulation and toxicological risk analysis. Nonetheless, accurate estimation of the drift and diffusion functions remains essential, and computational demands increase with the number of simulated paths.

Close

Figure 5:

The Milstein sample paths showing multiple stochastic trajectories of logKp versus time for a representative compound.

Export to PPT

Results of the Heston model [Figure 6] demonstrated improved handling of uncertainty, with predictive intervals aligning more closely to observed experimental spread. Evaluation metrics of the model included the mean absolute error (MAE), mean squared error (MSE), and the coefficient of determination (R²). The Heston model offers a promising pathway to incorporate stochastic dynamics into permeability predictions, bridging mathematics and percutaneous penetration. Optimal parameters (e.g., κ_y≈1.2, κv≈0.8, θ_v≈0.2) yielded predictive distributions consistent with observed scatter. The stochastic variance term captured a wider experimental spread in permeability. The Heston model achieved a mean standard error (MSE) of 0.66 and an R² of 0.72. Compared with Ridge regression, the Heston framework improved both accuracy and uncertainty calibration. Predictive intervals encompassed ~92% of observed test values, consistent with nominal 90% coverage. The Heston stochastic volatility model introduces a dual dynamic: mean reversion toward molecular-descriptor-based permeability predictions, and variance evolution capturing unexplained variability. This framework mirrors biological realities: skin barrier variability, localized differences in stratum corneum thickness, and fluctuating hydration all introduce heteroskedasticity.

Close

Figure 6:

Observed versus predicted logKp values through the Heston model.

Export to PPT

Figure 7 compares the predicted versus true Log Kp values obtained using a Gradient Boosting Regression (GBR) model. The points are distributed along a diagonal trend, suggesting a generally good correlation between model predictions and experimental (true) values. The GBR model successfully captures the nonlinear relationships among molecular descriptors. Minor deviations from the diagonal line imply some residual error for compounds with extreme Log Kp values (i.e., highly lipophilic or hydrophilic compounds). Overall, the model exhibits robust predictive performance and generalization capability, typical of ensemble tree-based methods that reduce bias through boosting successive weak learners. The GBR model effectively balances variance and bias through its boosting iterations, making it suitable for QSPR modeling of permeability. Figure 8 presents a similar scatter plot comparing actual (experimental) and predicted Log Kp values using an SVR model. The SVR algorithm employs an RBF kernel that aims to fit the data within a defined error margin (ε-insensitive tube). Compared to GBR, SVR appears to slightly underpredict at higher Log Kp values, which may result from kernel smoothing effects or suboptimal hyperparameter tuning. Figure 9 compares observed (experimental) and predicted Log Kp values obtained through non-LR using the LMA–a hybrid approach combining the gradient descent and Gauss– Newton methods. The red dashed line represents the ideal correlation (y = x). The tighter clustering of points near the central line suggests good predictive precision for mid-range permeability compounds but slightly reduced sensitivity at the extremes. Unlike machine-learning regressors, the LMA-based nonlinear fit emphasizes parametric optimization, typically using sigmoidal or exponential models derived from Fickian diffusion kinetics. The consistency across a wide range of Log Kp values suggests that the model captures underlying physicochemical phenomena (e.g., molecular size, lipophilicity, and polarity) influencing permeability.

Close

Figure 7:

LogKp predictions using the gradient boosting regression.

Export to PPT

Close

Figure 8:

Actual versus predicted LogKp values using the support vector regression.

Export to PPT

Close

Figure 9:

Observed versus predicted logKp values through the Levenberg–Marquardt algorithm.

Export to PPT

DISCUSSION

Permeability models published in the literature emphasize deterministic QSPR.^[8,88] Stochastic modeling aligns with emerging views that variability is not just noise but an intrinsic property of the drug absorption processes. We also used machine learning in our study. Machine learning seeks to employ algorithms to identify patterns within data samples, therefore creating models with generalization abilities that can “learn by analogy.”^[89] A different method has recently been utilized for transdermal data. Hamadeh et al. employed a systematic methodology to accurately estimate the skin permeability of a topically applied chemical, utilizing descriptors of the permeant, vehicle, and skin conditions, alongside mechanistic insights into permeability across a finely discretized dermal strata.^[90] The approach focuses on a direct computation of permeability derived from the mechanical characteristics of the permeation context.^[90] The computational method can be seamlessly incorporated into algorithms for parameter discovery, optimization, and sensitivity analysis.^[90]

Accurate forecasting of skin permeability is crucial for efficient transdermal drug delivery (TDD).^[91] Abdallah et al., addressed this critical necessity to enhance TDD.^[91] A dataset of 441 records for 140 substances with diverse LogKp values was examined.^[91] The descriptor computation yielded 145 relevant descriptors.^[91] Regression analysis was performed utilizing machine learning models, including RF, CatBoost, Multiple LR, XGBoost, LightGBM (LGBM), and Artificial Neural Networks (ANN).^[91] XGBoost, LGBM, and gradient boost models demonstrated enhanced predictive accuracy relative to alternative models.^[91] In the models, the key attributes with significant impact on skin permeability were hydrogen bond donors, hydrophobicity, topological polar surface area, and hydrogen bond acceptors.^[91]

Baba et al. executed several compelling experiments, compiling an extensive dataset (322 entries) of in vitro human skin permeability coefficients.^[92] From a chemical standpoint, this was a diverse set of penetrants. Subsequently, the authors created a dataset of several descriptors and analyzed the compounds through Gaussian process regression (GPR) and nonlinear SVR.^[92] The SVR model exhibited slight superiority over the GPR model, with root mean square error, squared correlation coefficient, and MAE values of 0.29, 0.94, and 0.21, respectively.^[92] Zeng et al. developed a three-descriptor quantitative structure–activity/toxicity relationship model to evaluate the skin permeability of 274 compounds.^[93] The authors used an SV machine (SVM) with a genetic algorithm (GA) to analyze the dataset.^[93] The investigators evaluated a training set of 139 chemicals and reported that the SVM model had a root mean square error of 0.253 and a coefficient of determination of 0.946.^[93] A root mean square error of 0.302 and an R² of 0.872 were obtained for a test set of 135 compounds.^[93]

Rezaei et al. developed a 3D-QSPR model to evaluate the influence of the physicochemical characteristics of 211 compounds on skin permeability.^[94] The Kennard-Stone method was used to divide the dataset into a test set of 52 molecules and a training set of 159 molecules.^[94] GA, fractional factorial design, and successive projection algorithm (SPA) were utilized to ascertain the most important three-dimensional molecular structure.^[94] The indices chosen through different feature selection techniques were correlated with skin permeability parameters utilizing SVM and partial least squares methods.^[94] The SPA-SVM model produced a correlation coefficient of 0.96, a Q2 of 0.73, and an R² of 0.76.^[94] The hydrogen bonding acceptor and donor properties of the analyzed compounds may impact their skin transport.^[94] Interestingly, it was shown that permeability increased with increased hydrophobicity and diminished with higher MW.^[94] Hydrophobicity in the target molecule, as well as the configuration and location, exerted a significant impact on skin permeability.^[94]

A recent study sought to assess the effectiveness of nonlinear and linear models in forecasting the skin penetration of 11 compounds.^[95] The authors utilized machine learning based on the Gaussian method as well as QSPR to forecast the permeability coefficients of the compounds.^[95] Transdermal absorption was evaluated using Franz cells and a standardized protocol that incorporated high-performance liquid chromatography.^[95] A statistical analysis comparing expected and experimentally derived values was performed using MSE and the Pearson sample correlation coefficient.^[95] The study’s findings revealed that the models poorly anticipated permeation, with certain instances deviating by two or three orders of magnitude from the experimentally obtained values.^[95] Nonetheless, this assortment of compounds allowed the models to precisely rank the permeants.^[95] The scientists concluded that although the models are unsuitable for accurate permeation predictions, they can successfully establish a rank order of permeation, facilitating the selection of candidate compounds for in vitro screening.^[95] It is crucial to recognize that these predictions must take into account the true relative potencies of pharmaceutical alternatives.^[95]

Lindh et al. utilized an aggregated conformal prediction (ACP) framework to predict the skin transport (log Kp) of compounds.^[96] ACP offers a prediction interval at a set confidence level for each molecule, enhancing informed decision-making, such as suggesting the next product for synthesis.^[96] The models were developed using SVM and RF techniques, based on empirically derived skin transport data from 211 different compounds.^[96] The resulting models exhibited comparable predictive parameters compared to previously published models, while additionally offering a reliable, individually defined prediction range for each chemical, rather than a single predicted value.^[96] The approach employs calculated indices and can swiftly predict the skin penetration rate of new chemicals.^[96]

Gradient boosting (GB) is a potent technique used in machine learning for optimization.^[26] Leo Breiman created the technique as a “arcing method” to reduce variability in classification.^[26] Jerome Friedman also employed classification and regression.^[26] GBR develops models incrementally, with each new model addressing the deficiencies of its forerunners.^[97] This approach often yields improved predictive precision. Gradient boosting regression (GBR) amalgamates the strengths of multiple weak learners, especially DTs, to alleviate overfitting and improve generalization.^[97] In machine learning, GBR adeptly models complex, non-linear relationships between targets and features, making it more versatile and effective for diverse datasets.^[97] The method includes factors such as tree complexity and learning rate, which help mitigate overfitting and ensure effective generalization to new data.^[97] Gradient boosting regression (GBR) can be used to handle various data types, including numerical and categorical variables, making it a suitable choice for numerous applications.^[97]

SVR is an efficient machine learning method that uses SVs to determine a hyperplane that best accommodates the data while minimizing the margin.^[98] SVR can be used to investigate nonlinear data and has diminished susceptibility to outliers relative to LR.^[98] SVR prioritizes the overarching trend of the data, disregarding specific outliers that significantly deviate, hence improving predictions.^[98] It uses several kernel functions such as sigmoid, polynomial, linear, and RB.^[98] This adaptability enables the modeling of diverse data patterns and can improve forecast accuracy.^[98] The LMA is a non-LR technique utilized in ANN to tackle various problems.^[99] The LMA is a conventional optimization method that combines the Gauss–Newton approach with the gradient descent methodology.^[99] The calculation of the Jacobian matrix is essential in the LMA, affecting its convergence and performance.^[99]

Three regression models were utilized to evaluate the relationship between permeability coefficient (target) and different features – MW, partition coefficient, Balaban distance connection index, Harary index, and Forgotten index. We applied machine learning methodologies to a revised Flynn dataset. We introduced the Forgotten index as an independent variable. After applying GBR, the R² was 0.84, and the mean square error (MSE) was 0.45. SVR with RBF kernel is a non-linear model that maps input features into a higher-dimensional space through a Gaussian kernel. This allows the model to discern complex, non-linear relationships between targets and features. The RBF kernel is commonly utilized in SVR because it may model many relationships without requiring the specification of the non-linearity’s structure. In the implementation of SVR, the R² and mean squared error (MSE) were measured at 0.78 and 0.62, respectively. The SVR model was created to predict the logarithm of skin permeability (logkp) using six aforementioned features. The SVR model utilized an RBF kernel, suitable for nonlinear relationships. The R² and MSE values were 0.58 and 0.94, respectively, following the use of the LMA. This method is an optimization method designed for non-linear least squares problems. The efficient identification of a model from a large dataset is a function of parameter optimization.^[36] Gradient-based optimization methods are commonly utilized for nonlinear systems because of their high computational efficiency and strong reliability.^[36]

The EM method enhances the predictive modeling of skin permeability by embedding random variability directly into the forecasting process. The presented figures – parity plot, residuals versus prediction, and predictive uncertainty distribution – collectively demonstrate that EM mean predictions align closely with experimental log Kp values, that residuals are largely unbiased, and that predictive uncertainty is well characterized. This approach provides pharmaceutical scientists and biomedical engineers with both accurate point estimates and rigorous quantification of uncertainty, thereby supporting more informed decisions in the design of TDD systems. We also used the Milstein approach. By incorporating the MS into QSPR modeling of skin permeability, researchers can move beyond point estimates to full probabilistic forecasts. This approach enhances the robustness of permeability predictions and provides a quantitative framework for uncertainty, aligning modeling practices with the inherent stochasticity of skin transport processes. The Heston method provides a rigorous framework for probabilistic prediction of skin permeability. Theoretical predictions regarding the cutaneous transport of chemicals are crucial for the development of transdermal therapeutic systems as well as the risk assessment of dangerous substances. In this project, we have used methods from stochastic calculus and machine learning.

The accurate prediction of skin permeability coefficients using SDE methods (EM algorithm, MS, and Heston model) and machine learning approaches (LMA, gradient boost regression, and SVR) holds significant potential for optimizing TDD systems across various therapeutic indications. In the context of chronic pain management, for instance, these predictive models can guide the rational design of transdermal patches by enabling formulators to screen candidate analgesic compounds based on their predicted permeation characteristics, thereby identifying molecules with optimal physicochemical properties for passive diffusion through the stratum corneum. The stochastic methods, particularly the EM algorithm with its superior predictive performance (R² = 0.81), can account for the inherent biological variability in skin permeability arising from differences in skin hydration, lipid composition, and regional anatomical variations among patients, thus providing more realistic permeability estimates that inform dosing regimen design. Meanwhile, machine learning approaches such as gradient boosted regression, which demonstrated the highest predictive accuracy (R² = 0.84), can rapidly evaluate extensive compound libraries during early-stage formulation development, enabling pharmaceutical scientists to prioritize lead candidates and optimize excipient selection to enhance drug penetration. These predictive tools facilitate the determination of appropriate drug loading concentrations, patch surface areas, and release kinetics required to achieve therapeutic plasma levels while minimizing local irritation and systemic toxicity. Furthermore, by correlating molecular descriptors such as MW, partition coefficient (log Kow), and topological indices with permeability outcomes, formulators can strategically modify drug molecules through prodrug approaches or select appropriate permeation enhancers that temporarily and reversibly alter skin barrier function. The integration of these computational methods into pharmaceutical development workflows ultimately accelerates the formulation optimization process, reduces the need for extensive in vitro and in vivo permeation studies, and supports the development of patient-friendly transdermal delivery systems for indications ranging from hormone replacement therapy and smoking cessation to the management of cardiovascular diseases and neurological disorders.

CONCLUSION

In this research project, the EM algorithm, the MS, and the Heston stochastic model were used to predict the permeability coefficient. Following the use of the EM method, the MSE was 0.5259, and the R² was 0.81. When the Milstein method was used, the MSE was 0.634, and the R² was 0.74. When we employed the Heston model, the mean standard error (MSE) was 0.66, and the R² was 0.72. Three other regression techniques (gradient boosting regression, SVR, and the LMA) were also employed to analyze skin permeation data. When gradient boost regression was utilized, the MSE was 0.45 while the R-squared was 0.84. Following the use of SVR, the mean squared error (MSE) was 0.62, and R² was 0.78. The application of the LMA resulted in an MSE of 0.94 and an R² of 0.58. Optimization of machine learning algorithms augments predictive accuracy and diminishes the time, resources, and labor expended on skin transport studies. The integration of SDE methods (particularly EM) with machine learning approaches (especially Gradient Boost Regression) provides a comprehensive framework for skin permeability prediction that balances predictive accuracy with the ability to model inherent biological variability, ultimately supporting rational TDD system design.