Machine Learning Models to Accelerate the Design of Polymeric Long-Acting Injectables

doi:10.21203/rs.3.rs-1623191/v1

Long-acting injectables are considered one of the most promising therapeutic strategies for the treatment of chronic diseases as they can afford improved therapeutic efficacy, safety, and patient compliance. The use of polymer materials in such a drug formulation strategy can offer unparalleled diversity owing to the ability to synthesize materials with a wide range of properties. However, the interplay between multiple parameters, including the physicochemical properties of the drug and polymer, make it very difficult to intuitively predict the performance of these systems. This necessitates the development and characterization of a wide array of formulation candidates through extensive and time-consuming in vitro experimentation. Machine learning is enabling leap-step advances in a number of fields including drug discovery and materials science. The current study takes a critical step towards data driven drug formulation development with an emphasis on long-acting injectables. A series of machine learning algorithms were trained and refined for accurate prediction of experimental drug release profiles. Analysis of the best performing model uncovered the properties of the drug and polymer that were identified to be key determinants of drug release. This information can be used to identify promising drug-polymer combinations that result in long-acting injectables with specific drug release behaviour. Importantly, implementation of this data driven approach has the potential to reduce the time, cost and resources associated with formulation development. Datasets and relevant codes used to train the machine learning models have been made openly available to encourage implementation in future drug formulation efforts.

long-acting injectables

machine learning

drug delivery

materials science

Long-acting injectables (LAI) are a class of advanced drug delivery systems that are designed to release their cargo over extended periods of time in order to achieve a prolonged therapeutic effect. LAIs that are amenable to parenteral administration can confer several advantages over conventional drug formulations, including increased patient compliance and bioavailability of drug.¹ Moreover, LAIs can be engineered to provide either local (e.g., Zilretta®) or systemic (e.g., Lupron Depot®) drug exposure over a prolonged period, making them ideal formulation strategies for the treatment of chronic diseases.² The unparalleled chemical and physical diversity afforded by polymers, makes polymer-based LAIs a particularly apt version of this drug delivery strategy. These systems can be engineered to entrap drug within a polymer matrix with release occurring via various mechanisms including erosion, diffusion, simultaneous erosion and diffusion, etc.³ In addition to achieving sustained or controlled drug release, the encapsulation of drug into these polymeric matrices can often protect therapeutic cargo.⁴ To date, various polymer-based LAI technologies administered via the intramuscular,⁵ subcutaneous,⁶ and intra-articular⁷ routes have received regulatory approval (Scheme 1A).

Despite the advantages associated with polymeric LAIs, their translation from bench to bedside remains non-trivial. In the past two decades, only about 30 polymeric LAI products have received regulatory approval, in contrast to the thousands of conventional oral formulations approved in the same period.^3,8 Several pivotal challenges limit the clinical translation of polymeric LAIs. Firstly, there are few biodegradable polymer materials that are generally recognized as safe (GRAS) for parenteral administration. To date, the polymeric LAIs which have received clinical approval are largely based on a single polymer: poly(lactide-co-glycolide) (PLGA).^2,3,9 The use of a material, such as PLGA, with an established safety profile may accelerate the regulatory approval of new LAI formulations. However, it is well recognized that polymer-drug compatibility significantly influences the performance of a formulation including drug loading capacity, drug release, and stability.¹⁰ Given that each drug has its own unique physicochemical properties it is unlikely that any one polymer material is ideally suited for the formulation of all drugs. Thus, reliance on the relatively small subset of polymeric materials that have GRAS status likely restricts our ability to develop polymer-based LAIs for many classes of drugs. Moreover, for a given polymer material there are a wide range of variables that must be optimized during LAI preparation. Moreover, changes to any of these variables has the potential to impact formulation performance, and the net effect of such alterations cannot be known a priori. Often, initially promising formulations can fail at various stages during development due to unwanted drug release rates, making their re-formulation and re-evaluation necessary (Scheme 1B). This trial-and-error based approach represents a significant bottleneck in LAI development.¹¹

To date, several strategies have been investigated to inform decision making and expedite the drug formulation development process. For instance, mathematical models have been used to describe and greatly enhance our understanding of drug release mechanisms.¹² However, the application of these empirical models is limited to post hoc analysis of LAI in vitro drug release profiles to determine mechanistic principles, and they do not offer information on in vitro drug release from LAIs a priori. More recently, molecular dynamics simulations have been investigated.¹³ These techniques have been useful in quantifying links between drug release rates and formulation parameters (including particle size and drug loading levels).¹⁴ While the development of these techniques is an active area of research, molecular scale simulations of entire drug delivery systems are currently computationally intensive. While these approaches can be used to confer useful information on potential LAI systems, they cannot currently be used in place of experimental drug release assays.¹⁴ Several studies have also investigated machine learning (ML) approaches.¹¹

The application of ML in the pharmaceutical sciences is generally limited by a lack of available open source datasets to train models.¹¹ Past efforts to predict in vitro drug release from LAIs using ML have exclusively considered neural network (NN) based models, and have examined narrow application domains. For example, Szlęk et al. used NNs to predict drug release of proteins and peptides from PLGA-based MPs using a dataset of 68 PLGA formulations extracted from the literature.¹⁵ Small molecule studies have also focused on NNs to predict drug release and these have generally been limited to less than 20 data instances.¹⁶ The use of NNs for supervised learning tasks in the low-data regime may be associated with increased risk of overfitting compared to alternative ML algorithms, such as tree-based models or Gaussian processes, which are better suited for sparse data problems.

The current study pursued the investigation and refinement of ML models for accurate prediction of fractional drug release from polymeric LAIs (Scheme 1C). To achieve this, we trained and evaluated a series of eleven different ML algorithms, including multiple linear regression (MLR), MLR with least absolute shrinkage and selection operator regularization (LASSO), partial least squares (PLS), decision tree (DT), random forest (RF), light gradient boosting machine (LightGBM), extreme gradient boosting (XGBoost), natural gradient boosting (NGBoost), support vector regressor (SVR), k-nearest neighbors (k-NN), and neural networks (NN). The refined RF model was found to have the best performance: predicting fractional drug release with a high degree of accuracy and furnishing these predictions with estimates of uncertainty.

Dataset development

Dataset collection

Two datasets, (i) a training dataset, and (ii) an external validation dataset, were constructed from previously published studies. Each dataset comprised a range of different small molecule drugs as well as both MPs and cylinder-shaped polymer LAI systems. The LAIs were formed from several commercially available polymers, such as PLGA, polylactic acid (PLA), and polycaprolactone (PCL), of various molecular weights and lactide to glycolide ratios, for PLGA. The training dataset used for ML model screening and development contained 102 drug release profiles for 22 unique drug-polymer combinations. In total this comprised a total of 2638 individual fractional release measurements. This data was either generated in-house and previously published by our research group or collected from previously published studies. Data obtained from external sources was collected using Web of Science and the keyword combinations of “polymeric microparticle” and “drug delivery”. The second dataset (i.e., the external validation dataset) was used to conduct a pseudo-prospective study on the trained ML models with the optimal hyperparameter configurations identified using the training dataset. The external validation dataset was exclusively collected from an additional Web of Science search and contained 79 unseen release profiles from 21 drug-polymer combinations for a total of 1016 data entries (Table 1). A visual representation of both datasets with respect to their target values is shown in Figure S1.

Table 1

Summary of drugs, polymers, and LAI properties for formulations included in the external validation dataset. *LAI size refers to the average diameter for spherical particles and average width (W) and height (H) dimensions for cylinders.
Drug	Polymer molecular weight (average)	Drug loading capacity	LAI size* (average)	LAI shape
Dexamethasone (DEX)	PEA (49 kD) PLGA (28 kD; L/G = 50/50) PLGA (13 kD; L/G = 50/50)	15% 8% 8%	14 µm 20 µm 20 µm	Spherical
Gefitinib (GEF)	PLGA (44 kD; L/G = 50/50)	3–8%	5–130 µm	Spherical
Cannabidiol (CBD)	PLGA (12 kD; L/G = 50/50) PLGA (51 kD; L/G = 50/50) PCL (42.5 kD)	8% 5% 9%	24 µm 24 µm 50 µm	Spherical
Tetrahydrocannabinol (THC)	PCL (42.5 kD)	8%	50 µm	Spherical
Caffeine (CAF)	PLGA (46 kD; L/G = 50/50)	6%	62–94 µm	Spherical
Lidocaine (LDC)	PLGA (12 kD; L/G = 50/50) PLGA (48 kD; L/G = 50/50) PLGA (10 kD; L/G = 70/30)	1–18% 1–2% 5–42%	3–9 µm 4 µm 5–100 µm	Spherical
Etoricoxib (ETC)	PCL (10 kD)	3%	16 µm	Spherical
Quercetin (QRC)	PCL (14 kD)	3–4%	61–100 µm	Spherical
Paclitaxel (PTX)	PLA (106 kD) PLGA (57.5 kD; L/G = 50/50) PLGA (10.5 kD; L/G = 75/25) PVL-co-PAVL (15 kD) PVL-co-PAVL (32 kD) PVL-co-PAVL (39 kD)	2% 1–2% 1–2% 10% 10–20% 10–20%	0.6 µm 0.4–0.6 µm 0.3–0.7 µm 3 mm×2 mm (W×H) 3 mm×2 mm (W×H) 3 mm×2 mm (W×H)	Spherical Spherical Spherical Cylinder Cylinder Cylinder
Triamcinolone acetonide (TAA)	PLGA (18 kD; L/G = 50/50) PLGA (54 kD; L/G = 50/50) PVL-co-PAVL (15 kD) PVL-co-PAVL (32 kD) PVL-co-PAVL (39 kD)	5% 5% 10% 10–20% 10%	73 µm 71 µm 3 mm×2 mm (W×H) 3 mm×2 mm (W×H) 3 mm×2 mm (W×H)	Spherical Spherical Cylinder Cylinder Cylinder
Triamcinolone hexacetonide (TAH)	PVL-co-PAVL (15 kD) PVL-co-PAVL (39 kD)	10% 10%	3 mm×2 mm (W×H)	Cylinder
Curcumin (CCM)	PVL-co-PAVL (15 kD)	10%	3 mm×2 mm (W×H)	Cylinder
Acetaminophen (ACE)	PVL-co-PAVL (15 kD) PVL-co-PAVL (39 kD)	10%	3 mm×2 mm (W×H)	Cylinder
Letrozole (LTZ)	PLGA (57.7 kD; L/G = 65/35)	12–20%	0.6–1 µm	Spherical
Diprophylline (DPP)	PLGA (46 kD; L/G = 50/50)	5–7%	63–296 µm	Spherical
Ibuprofen (IBP)	PLGA (46 kD; L/G = 50/50)	14–20%	91–383 µm	Spherical

Initial feature selection

Fourteen molecular and physicochemical descriptors based on domain knowledge were initially selected as input features to describe the LAI formulations to the various ML models. These included the drug’s molecular weight (Drug_MW), topical polar surface area (Drug_TSPA), number of heteroatoms (Drug_NHA), melting temperature (Drug_Tm), acid dissociation constant (Drug_pK_a), and calculated partition coefficient (Drug_LogP), the polymer’s molecular weight (Polymer_MW), molar cross-linking ratio (CL_Ratio), and lactide-to-glycolide ratio (LA/GA), the drug loading capacity of the LAI expressed as a fraction (DLC), the surface area-to-volume ratio of the LAI (SA-V), the percent of surfactant dissolved in the experimental release media (SE), and drug release time points in days (Time) for each fractional drug release measurement. With the exception of SA-V, the 14 input features were either extracted from original publications or calculated using the RDkit package. SA-V was constructed and implemented for this study as it confers information that is related to the size and shape of the LAI system. This enables the inclusion of both spherical and cylindrical shaped LAIs in one model.

Comparison of different ML algorithms

Model training and cross-validation

Each ML model was subject to a hyperparameter optimization procedure using leave-one-group-out (LOGO) cross-validation, in which fractional release values for each drug-polymer combination in the training dataset constituted a single group. The range of model hyperparameters considered, as well as the final selected values are shown in Table S1-7. Model hyperparameters are tuned using a random grid search, where the objective function was the average model performance across each held-out group. Model performance was assessed using the mean absolute error (MAE) metric, which is the average absolute difference between predicted and experimental fractional release values. MAE values achieved using final hyperparameter settings are displayed in Table 2(a). We found that tree-based models were on average more accurate (MAE < 0.2) than the linear, instance-based, and deep learning models investigated here. The RF model had the highest overall prediction accuracy, with MAE = 0.164

Pseudo-prospective study

To evaluate the generalizability of the trained ML models a pseudo-prospective study was conducted using the external validation dataset. The results of the pseudo-prospective studies are summarized in Fig. 1. Overall, the predictive performance of these models on the external validation dataset largely reflected their respective accuracies on the training dataset. Specifically, the tree-based ML models were more accurate than the other model classes, with the RF model resulting in the lowest MAE value (Table 2(b)). Apart from the RF model, the performance of the other models generally regressed in the pseudo-prospective study. For instance, the MAE of the DT model increased from 0.182 on the training dataset to 0.204 on the external validation dataset, while the MAE of the SVR model increased from 0.234 to 0.301. Interestingly, the RF model was closely matched in performance by another, relatively new, tree-based model; NGBoost. The NGBoost algorithm was first reported in 2019, and it employs natural gradients to afford probabilistic predictions.¹⁷ While, the performance of the RF model was not significantly better than the NGBoost model (p = 0.08), the RF model was selected for further development and analysis in this study.

Table 2. Mean absolute error (MAE) values obtained for the different ML algorithms following prediction of fractional drug release for (a) the training dataset, and (b) the external validation dataset (pseudo-prospective study). MAE values in the table are color coded as follows: MAE £ 0.2 highlighted in green; MAE values > 0.2 and £ 0.3 are highlighted in yellow; and MAE values > 0.3 are highlighted in red. The standard error of the mean (σM) (Eq. 3) and median absolute error (Median AE) results are also shown for the pseudo-prospective study.

(a) MAE values obtained on training dataset following LOGO cross validation (n = 22 folds).
	RF	NGBoost	XGBoost	LightGBM	DT	PLS	k-NN	NN	LASSO	SVR	MLR
MAE	0.164	0.170	0.181	0.167	0.182	0.311	0.210	0.218	0.330	0.234	0.359
(b) MAE values obtained during the pseudo-prospective study (n = 1017 samples)
	RF	NGBoost	XGBoost	LightGBM	DT	PLS	k-NN	NN	LASSO	SVR	MLR
MAE	0.164	0.172	0.186	0.191	0.204	0.210	0.226	0.229	0.252	0.301	0.383
σM	0.004	0.004	0.005	0.005	0.005	0.005	0.005	0.006	0.005	0.009	0.010
Median AE	0.126	0.133	0.158	0.156	0.164	0.179	0.199	0.172	0.249	0.233	0.311

RF model refinement

Feature engineering

The RF model identified as the most accurate at predicting fractional drug release during the pseudo-prospective study consisted of the 14 input features mentioned above. Input feature sets for ML models should be comprehensive and encode all relevant physical information to make accurate predictions, adding too many, possibly irrelevant features can increase the risk of introducing spurious correlations between individual features and observed release profiles which can eventually degrade the model’s generalizability. Thus, the 14-feature RF model was assessed, refined, and optimized through agglomerative hierarchical clustering analysis to identify potentially redundant input features. The 14 input features were arranged into a hierarchy of clusters using the farthest neighbor clustering algorithm (Fig. 2a). The performance of the RF model was then assessed following the removal of select clusters of features based on their linkage distances (Fig. 2b). The optimal version of the RF model was determined to be the model that could achieve the lowest MAE while minimizing the number of features. This was identified to be a 12-feature version of the RF model that did not include Drug_NHA or Drug_TPSA. It was found that the removal of any additional feature clusters resulted in poorer performance. Thus, it appears that the strong correlation observed between Drug_NHA, Drug_TPSA, and Drug_MW resulted from some redundancy in the 14-feature RF model (Fig. 2).

Confidence intervals for the RF model

A method previously developed by Wager et al. and Efron et al. that employs the jackknife and infinitesimal jackknife resampling for variance estimation was used to generate uncertainty estimates for the 12 feature RF model predictions.^18–20 Fig. 3 and Figure S5-6, show the predicted drug release profiles and associated uncertainty estimates in comparison to the experimental fractional drug release profiles. Overall, the predicted and experimental drug release profiles showed a high degree of similarity. In many cases, where the fractional drug release predictions did not match well with the experimental release values, the variance extracted from the RF model could account for these differences. Relevant examples include PTX-PVL-co-PAVL, and IBP-PLGA predictions shown in Fig. 3, as well several examples in the supplementary information (i.e., CAF_PLGA_index10, ETC-PCL_index15, and LTZ-PLGA_index38 (Figure S5-6)). ML models that return uncertainty estimates can afford users more confidence in predictions, and this can be useful in applications (such as active learning) to identify the most informative data to collect in order to refine predictions.²¹

RF model analysis

RF feature importance analysis with SHAP

The effects of the various input features were analyzed using SHAP analysis on the 12 feature RF model. Figure 4A and Fig. 4B shows the global contribution for each feature in terms of SHAP value and absolute SHAP value, respectively. In each panel of Fig. 4, the input features are arranged from top to bottom in order of decreasing impact on prediction. For instance, Time is listed in the top row, indicating that it had the greatest impact on the fractional drug release predictions. In addition, each row of the swarm plot in Fig. 4A also depicts an individual dot plot, with each dot representing a separate instance in the external validation dataset. Blue dots indicate lower values for that feature in the dataset, while the red dots represent higher values. For example, blue dots in the Time row represent low values for drug release time points, while pink dots in the Polymer_MW row represent instances for LAIs composed of polymers with higher molecular weights in the external validation dataset. The distance between an individual dot and the grey vertical central line reflects the extent to which that fractional drug release value deviates from the mean value of the external validation dataset. For example, lower Time values (i.e., blue dots in the time row, Fig. 4A) resulted in a strong negative contribution to fractional drug release, while higher Time values (i.e., red dots in the time row, Fig. 4A) resulted in a strong positive contribution. After Time, Polymer_MW and LA/GA are the next most influential features. Interestingly, both are descriptors of the polymers used to construct the LAIs which indicates that, of the components of the LAI system, the RF model has recognized that the polymer properties have the most significant influence on fractional drug release. The remaining physicochemical and molecular descriptors appear to be largely irrelevant for the RF model. However, it should be noted that Fig. 4A (and Fig. 4B) only shows the effect of each individual feature and does not account for potential synergy between input features. To better visualize the collective effect of all 12 input features the decision path function from the SHAP library was utilized. Figure 4C shows the decision path taken by each prediction in the external validation dataset. This type of analysis allows individual data instances to be followed and interpreted.

Decomposition of RF model predictions with SHAP

Figure 5 shows a select example of an LAI formulation taken from the external validation dataset to illustrate how the fractional drug release predictions are generated by the 12 feature RF model. LTZ-PLGA Experimental_index39 is a spherical MP (1 mm) of letrozole loaded into PLGA (L/G = 63:35, Mw = 57.7 kD) at a DLC of 17% (wt%). Figure 5A shows a visual comparison of the experimental and predicted fractional drug release values for LTZ-PLGA Experimental_index39, and highlights three arbitrary prediction instances (i) Time = ~ 2 days, predicted fractional drug release (i.e., f(x)) = 0.22; (ii) Time = ~ 8 days, f(x) = 0.39; and (iii) Time = ~ 20 days, f(x) = 0.65. Figure 5B shows the decision path for each of the predicted fractional drug release values. Each prediction begins at the base value for the model and migrates through the RF based on contributions of each input feature. Figure 5C(i-iii) shows a breakdown of the exact contribution of each input feature on the prediction of fractional drug release for the three highlighted data instances, as force plots. In Fig. 5C, the exact values of the input features are given, and the relative contribution of each feature is shown by a pink (positive) or blue (negative) band on the force plot, with the width of this band representing the numerical contribution to the final model output. The contributions of each input feature are summed with the base value of the model to derive the model output. This is further illustrated in Figure S2-4.

Feature interactions with SHAP

To investigate the relationship between Time and the other input features we constructed SHAP dependence plots (Fig. 6). The values of Time are plotted along the x-axis, the calculated SHAP values for Time are plotted on the y-axis, and the data distributions for select interaction features including LA/GA, Polymer_MW, DLC, and SA-V are shown along a second y-axis where they are colored according to their relative value in the external validation dataset. Figure 6B(i) and Fig. 6B(iii) show a fairly uniform distribution of LA/GA values across time values, indicating that there is no particular LA/GA that is more likely to release drug at a faster or slower rate. Figure 6B(iii) shows that only higher DLC values (i.e., a lack of blue data points in the highlighted region) are present at higher timepoints. Indicating that the model may recognize that LAIs with higher drug loading levels afford drug release over longer periods of time. Figure 6B(ii) show a cluster of high molecular weight polymers at higher time values, as highlighted by the black dashed circle. Figure 6B(iv) shows a similar cluster of LAIs associated with low surface-area to volume ratios (i.e., larger sized systems) at higher time values (also highlighted by a black dashed circle). This indicates that in our external validation dataset, only high molecular weight polymers and LAI systems with a large size released drugs at later timepoints. In this way, the results of SHAP analysis can provide general guidance on the design of LAIs for specific drug release applications. For example, to design an LAI system that releases drug for durations of more than 40 days, the model suggests systems formed from high molecular weight polymers (i.e., Polymer_MW > 40 kD) that are larger in size (i.e., SA-V < 2000 mm^-1) with a relatively high drug loading capacity (i.e., DLC \(\ge\)0.1 (10 wt%)).

Cutting-edge ML technologies are now freely available to pharmaceutical and material scientists, and the results obtained here demonstrate the potential for ML to expedite the development of innovative drug delivery technologies. Among the strengths of modern ML techniques are their ability to provide uncertainty estimates reflecting model confidence, and to provide insights into how the models reach their predictions. Herein we demonstrate that ML models can be used to predict in vitro drug release from LAIs with a high degree of accuracy despite a minimal amount of training data (i.e., ~ 100 LAI formulations). The implementation of these models has the potential to reduce the time and cost associated with future LAI formulations. This is promising for the application of ML in drug formulation development, and other industrial sectors. In particular, where large, reliable, open access experimental datasets are not readily available. This proof-of-concept study, and its associated dataset, will aid in fostering the development of more advanced, tailored, and accurate ML approaches for the prediction of drug release profiles a priori. This set the stage for a new era of data-driven development in advanced drug delivery research.

Data collection

Two datasets were constructed to train and validate ML models for fractional drug release prediction; (i) training dataset, and (ii) external validation dataset. The training dataset used for ML model screening and development. It consisted of 2638 data entries, which described 102 drug release profiles for 22 drug polymer combinations. The data was either generated in-house and previously published by our lab or collected from Web of Science. The data obtained from Web of Science was collected using the keyword combinations of “polymeric microparticle” and “drug delivery”. From appropriate manuscripts, information related to the preparation, final composition, and release kinetics of drug-loaded LAI were collected. The latter was primarily extracted from figures of drug release profiles using the “GetData Graph Digitizer” application. The external validation dataset was exclusively collected from an additional Web of Science search using the same search terms. This dataset consisted of 1016 data entries, which described the release profiles collected for 79 unseen LAI formulations made up from 21 drug-polymer combinations.

The initially collected dataset was composed of a table of drug and polymers names, as well as formulation physicochemical properties, and fractional drug release values at various timepoints. In order to use this data to construct and train ML models it is necessary to describe various elements using computational recognizable descriptors which were collected from RDkit (i.e., using SMILES codes). The polymers and LAI formulations were described exclusively using information reported in the relevant published articles, these included; polymer molecular weight, lactide-to-glycolide ratio (for non-PLGA systems this was set as zero), molecular crosslinking ratio of polymers (for non-cross-linked systems this was set as zero), initial drug-to-polymer ratio, drug loading capacity (DLC), surface area-to-volume (SA-V) ratio for the LAI system, and the precent of surfactant present in the release media (SE; where no surfactant was present in the release media, this was set as zero).

Data splitting strategy for ML model training

In all cases drug-polymer combinations were grouped as such, to allow cross-validation against drug-polymer based splits. Data splitting was done using the group-shuffle-split (GSS) method combined the leave-one-group-out (LOGO) cross-validation strategy. The GGS and LOGO packages were all imported from the Scikit-learn library in Python.²²

ML model development

In total eleven ML algorithms were trained and investigated for this task. These included, multiple linear regression (MLR), MLR with least absolute shrinkage and selection operator regularization (LASSO), partial least squares (PLS), decision tree (DT), random forest (RF), light gradient boosting machine (LightGBM), extreme gradient boosting (XGBoost), natural gradient boosting (NGBoost), support vector regressor (SVR), k-nearest neighbors (k-NN), and neural network (NN) models. All of these ML models were built and evaluated in Python. NN models were built using the Keras package with the backend of TensorFlow,²³ LighGBM models were built using the LightGBM package, XGBoost models were built using the XGBoost package,²⁴ NGBoost models¹⁷ were built using the NGBoost package, and all other models were built using the Scikit learn library.²²

In all cases, prior to training any of the ML models, a data preprocessing step was conducted to standardize the data prior to input into the ML models. This was done using the standard scalar package available in the Scikit learn library.²² ML model hyperparameters were tuned using the randomized grid search package in Scikit learn,²² and the negative mean absolute error metric was employed. The hyperparameters screened, and numerical range for these values, varied depending on the ML model being training. This screening process is summarized in Table S1-7.

Pseudo-prospective study (ML model evaluation)

In order to better assess the predictive performance of all the trained ML models a pseudo-prospective study was conducted using the external validation dataset. This dataset (~ 1000 data samples) enabled a more quantitative evaluation of the prediction accuracy of each model. This was done by determining the absolute error (AE) for each prediction made by each model, as well as the mean absolute error (MAE) for various groups (i.e., drugs, polymer, and drug-polymer combinations) within the external validation dataset. The AE, MAE, and standard error of the mean () were determined using the equations shown below

\(absolute error \left(AE\right)= \left|{y}_{i}-{x}_{i}\right|\)

Eq. 1

\(mean absolute error \left(MAE\right)= \frac{\sum _{i=1}^{n}\left|{y}_{i}-{x}_{i}\right|}{n}\)

Eq. 2

\(standard error of the mean \left(\sigma M\right)= \frac{s}{\surd n}\)

Eq. 3

Where, \({y}_{i}\) is the predicted fraction drug release value; \({x}_{i}\) is the experimental fractional drug release value obtained from either the training or external validation datasets; \(n\) is the total number of data points; and \(s\) is the standard deviation.

Feature engineering

Agglomerative hierarchical clustering analysis was performed using the hierarchical clustering package from SciPy in Python ²⁵ to arrange the initial input features into a hierarchy of clusters using the farthest neighbor clustering algorithm. The performance of the optimal ML algorithm (RF) was then assessed following the removal of select clusters based on their linkage distance. The hyperparameter structure for the RF model, that was identified during training (Table S8), remained fixed and only the number of input features were varied.

Confidence intervals for random forest model

Confidence intervals values for the optimal RF model were determined based on variance estimates for bagging originally proposed by Efron et al; the jackknife and infinitesimal jackknife estimators.^19,20 In this study, the variance for the optimal random forest model were calculated using the forestci package in Python and used to portray the confidence intervals for the fractional drug release predictions.²⁶

Model interpretation

Shapley additive explanation (SHAP) analysis was conducted on the 12 feature RF model trained for the prediction of fractional drug release from LAIs. The effect of the various input features on the fractional drug release prediction for the external validation dataset was assessed using the TreeSHAP package and the force plot visualizations from the SHAP library in Python. ^27,28

Acknowledgements

NSERC Discovery grant (RGPIN-2016-04293) to C.A. F.H., M.A., and A.A.G. acknowledge support by the Defense Advanced Research Projects Agency under the Accelerated Molecular Discovery Program under Cooperative Agreement No. HR00111920027 dated August 1, 2019. M.A. is supported by a Postdoctoral Fellowship from the Vector Institute.

Data Availability

The code, dataset, and results that support the findings of this study are available at the Aspuru-Guzik Group’s GitHub page (https://github.com/aspuru-guzik-group/long-acting-injectables). A preprint of this manuscript is also available on ChemRxiv (https://doi.org/10.26434/chemrxiv-2021-mxrxw-v2).

Brigham, N. C., Ji, R.-R. & Becker, M. L. Degradable polymeric vehicles for postoperative pain management. Nat. Commun. 12, 1367 (2021).
Ghitman, J., Biru, E. I., Stan, R. & Iovu, H. Review of hybrid PLGA nanoparticles: Future of smart drug delivery and theranostics medicine. Mater. Des. 193, 108805 (2020).
O’Brien, M. N., Jiang, W., Wang, Y. & Loffredo, D. M. Challenges and opportunities in the development of complex generic long-acting injectable drug products. J. Controlled Release 336, 144–158 (2021).
Siepmann, J. & Siepmann, F. Microparticles Used as Drug Delivery Systems. in Smart Colloidal Materials (ed. Richtering, W.) 15–21 (Springer, 2006). doi:10.1007/3-540-32702-9_3.
Agnihotri, S. M. & Vavia, P. R. Pharmacokinetics of intramuscular microparticle depot of valdecoxib in an experimental model. Drug Dev. Ind. Pharm. 35, 1043–1047 (2009).
Shenoy, D. B., D’Souza, R. J., Tiwari, S. B. & Udupa, N. Potential applications of polymeric microsphere suspension as subcutaneous depot for insulin. Drug Dev. Ind. Pharm. 29, 555–563 (2003).
Kraus, V. B. et al. Synovial and systemic pharmacokinetics (PK) of triamcinolone acetonide (TA) following intra-articular (IA) injection of an extended-release microsphere-based formulation (FX006) or standard crystalline suspension in patients with knee osteoarthritis (OA). Osteoarthritis Cartilage 26, 34–42 (2018).
Park, K. et al. Injectable, long-acting PLGA formulations: Analyzing PLGA and understanding microparticle formation. J. Controlled Release 304, 125–134 (2019).
Hines, D. J. & Kaplan, D. L. Poly (lactic-co-glycolic acid) controlled release systems: experimental and modeling insights. Crit. Rev. Ther. Drug Carrier Syst. 30, 257–276 (2013).
Bao, Z. et al. Poly(δ-valerolactone-co-allyl-δ-valerolactone) cross-linked microparticles: Formulation, characterization and biocompatibility. J. Pharm. Sci. 110, 2771–2777 (2021).
Bannigan, P. et al. Machine learning directed drug formulation development. Adv. Drug Deliv. Rev. (2021) doi:10.1016/j.addr.2021.05.016.
Irurzun-Arana, I., Rackauckas, C., McDonald, T. O. & Trocóniz, I. F. Beyond Deterministic Models in Drug Discovery and Development. Trends Pharmacol. Sci. 41, 882–895 (2020).
Salo-Ahen, O. M. H. et al. Molecular Dynamics Simulations in Drug Discovery and Pharmaceutical Development. Processes 9, 71 (2021).
Casalini, T. Not only in silico drug discovery: Molecular modeling towards in silico drug delivery formulations. J. Controlled Release 332, 390–417 (2021).
Szlęk, J., Pacławski, A., Lau, R., Jachowicz, R. & Mendyk, A. Heuristic modeling of macromolecule release from PLGA microspheres. Int. J. Nanomedicine 8, 4601–4611 (2013).
Li, Y., Rauth, A. M. & Wu, X. Y. Prediction of kinetics of doxorubicin release from sulfopropyl dextran ion-exchange microspheres using artificial neural networks. Eur. J. Pharm. Sci. 24, 401–410 (2005).
Duan, T. et al. NGBoost: Natural Gradient Boosting for Probabilistic Prediction. ArXiv191003225 Cs Stat (2020).
Wager, S., Hastie, T. & Efron, B. Confidence intervals for random forests: The jackknife and the infinitesimal jackknife. J. Mach. Learn. Res. 15, 1625–1651 (2014).
Efron, B. Jackknife-After-Bootstrap Standard Errors and Influence Functions. J. R. Stat. Soc. Ser. B Methodol. 54, 83–127 (1992).
Efron, B. Estimation and Accuracy After Model Selection. J. Am. Stat. Assoc. 109, 991–1007 (2014).
Reker, D. & Schneider, G. Active-learning strategies in computer-assisted drug discovery. Drug Discov. Today 20, 458–465 (2015).
Pedregosa, F. et al. Scikit-learn: Machine Learning in Python. J. Mach. Learn. Res. 12, 2825–2830 (2011).
TensorFlow Developers. TensorFlow. (Zenodo, 2022). doi:10.5281/ZENODO.4724125.
Chen, T. & Guestrin, C. XGBoost: A Scalable Tree Boosting System. in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining 785–794 (ACM, 2016). doi:10.1145/2939672.2939785.
Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261–272 (2020).
Polimis, K., Rokem, A. & Hazelton, B. Confidence Intervals for Random Forests in Python. J. Open Source Softw. 2, (2017).
Lundberg, S. M. et al. From local explanations to global understanding with explainable AI for trees. Nat. Mach. Intell. 2, 56–67 (2020).
Lundberg, S. M. et al. Explainable machine-learning predictions for the prevention of hypoxaemia during surgery. Nat. Biomed. Eng. 2, 749–760 (2018).

Scheme 1 is available in the Supplementary Files section

There is NO Competing Interest.

floatimage1.png
Scheme 1. A) Scheme demonstrating selected routes of administration for FDA approved LAI formulations. B) Schematic depiction of a typical trial-and-error loop commonly employed during the development of LAIs, termed “Traditional LAI formulation development”. C) Summary of the workflow employed here to train various ML models to predict drug release from LAIs, termed “Data-driven LAI formulation development”.
SupplementaryInformation.docx
Supplementary Information File

Machine Learning Models to Accelerate the Design of Polymeric Long-Acting Injectables

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Results And Discussion

Dataset development

Dataset collection

Initial feature selection

Comparison of different ML algorithms

Model training and cross-validation

Pseudo-prospective study

RF model refinement

Feature engineering

Confidence intervals for the RF model

RF model analysis

RF feature importance analysis with SHAP

Decomposition of RF model predictions with SHAP

Feature interactions with SHAP

Conclusion

Methods

Declarations

References

Scheme

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1