Interaction of systemic drugs causing ocular toxicity with organic cation transporter: an artificial intelligence prediction

Abstract Chronic disease patients (cancer, arthritis, cardiovascular diseases) undergo long-term systemic drug treatment. Membrane transporters in ocular barriers could falsely recognize these drugs and allow their trafficking into the eye from systemic circulation. Hence, despite their pharmacological activity, these drugs accumulate and cause toxicity at the non-target site, such as the eye. Since around 40% of clinically used drugs are organic cation in nature, it is essential to understand the role of organic cation transporter (OCT1) in ocular barriers to facilitate the entry of systemic drugs into the eye. We applied machine learning techniques and computer simulation models (molecular dynamics and metadynamics) in the current study to predict the potential OCT1 substrates. Artificial intelligence models were developed using a training dataset of a known substrates and non-substrates of OCT1 and predicted the potential OCT1 substrates from various systemic drugs causing ocular toxicity. Computer simulation studies was performed by developing the OCT1 homology model. Molecular dynamic simulations equilibrated the docked protein-ligand complex. And metadynamics revealed the movement of substrates across the transporter with minimum free energy near the binding pocket. The machine learning model showed an accuracy of about 80% and predicted the potential substrates for OCT1 among systemic drugs causing ocular toxicity – not known earlier, such as cyclophosphamide, bupivacaine, bortezomib, sulphanilamide, tosufloxacin, topiramate, and many more. However, further invitro and invivo studies are required to confirm these predictions. Communicated by Ramaswamy H. Sarma

The benefits of the drugs to treat chronic diseases or the lack of other therapies often justify the associated risk of ocular toxicities; thereby, patients usually end up continuing their systemic medication (Brock et al., 2013;Vishnevskia-Dai et al., 2021).However, patients suffering with irreversible vision loss are left with no alternative but to stop their treatment.Also, such toxicities vary among patients, making it even more challenging for physicians to provide timely intervention (Shin et al., 2020;Yuan et al., 2019).Many studies indicate a regular ophthalmic examination of patients (management among experts from different specializations) consuming anti-neoplastic, central nervous system (CNS), cardiovascular, anti-arthritis drugs, and any other drug used for a chronic period (Ali et al., 2022).Therefore, it is crucial to understand the entry of systemic drugs into the eye during the drug development stage to avoid future toxicities.
Systemic drugs selectively bind and accumulate in ocular tissues such as the conjunctiva, cornea, lens, choroid, or retina despite the ocular barriers-causing ocular toxicity.We hypothesize that the transporters present in ocular barriers for transporting endogenous molecules such as nutrients, vitamins, and neurotransmitters could falsely recognize the systemic drugs and enable their entry into the eye (Kubo et al., 2014;Nirmal et al., 2010;Nirmal et al., 2013).In the last few decades, there has been increasing evidence of the role of membrane transporters in drug accumulation at the non-target site leading to toxicity, such as taurine transporter for vigabatrin uptake in the eye, and organic cation transporter (OCT) 2 for cisplatin uptake in the kidney (Filipski et al., 2009;Police et al., 2020).It has also been highlighted as a scientific session, 'Transporters and Toxicity,' at the International Transporter Consortium (ITC) Workshop IV in 2021 (Hafey et al., 2022).Various ocular tissues show the expression of organic cation transporter in the following order (high to low), retina (Zhang et al., 2008), iris-ciliary body, the cornea (Garrett et al., 2008), conjunctiva (Garrett et al., 2008), and lacrimal gland (Velpandian et al., 2012).It is also highly expressed in the liver, followed by the kidney and small intestine (Zhang et al., 2008).Around 40% of commonly prescribed drugs are organic cations at physiological pH, most acting as substrates or inhibitors for these transporters (Neuhoff et al., 2003).Therefore, we aim to identify the potential substrates of OCT1 among the systemic drugs causing ocular toxicity to understand their entry into the eye despite the ocular barriers.
Screening thousands of drugs for their interaction with OCT1 using conventional invitro and invivo models is complex.However, artificial intelligence methods such as machine learning and computer simulations can aid in better understanding drug-transporter interactions and drug toxicity with high reproducibility and reliability (Jain & Ecker, 2019;Jensen et al., 2021;Khuri & Deshmukh, 2018;Liu et al., 2016;Pu et al., 2019;Vamathevan et al., 2019).In silico predictions of drug toxicity have also been approved by The Organisation for Economic Co-operation and Development (OECD) and other regulatory boards (OECD, 2020(OECD, , 2021)).Earlier studies performed by other groups and us have predicted potential OCT1 substrates among Food and Drug Administration (FDA) approved drugs from 2014 to 2019 using the artificial intelligence method.(Baidya et al., 2020;Hendrickx et al., 2013;Jensen et al., 2021).However, the potential OCT1 substrates among the systemic drugs causing ocular toxicity were not screened earlier.Therefore, in the present study, we applied an advanced approach of computer simulations to understand the drug-OCT1 interactions.Artificial intelligence models were developed based on the structural and physicochemical properties as well as the biological activity (OCT1 uptake ratio) to predict the substrate or non-substrates of OCT1.Molecular dynamics (MD) simulations were used to evaluate drug-OCT1 binding stability and interactions at the atomic level, whereas metadynamics simulations were used to visualise drug movement across the transporter (OCT1).

Dataset preparation
To develop the machine learning model for the human OCT1 (hOCT1) training dataset was prepared based on Human Embryonic Kidney (HEK-293) in vitro experimental data to determine the uptake of drugs by hOCT1 (Hendrickx et al., 2013).The drugs with an uptake ratio greater than 1.3 were considered substrates for hOCT1.The training dataset comprised 110 substrates and 84 non-substrates, as used earlier by our group (Baidya et al., 2020).The screening dataset consisted of more than 600 systemic drugs causing ocular toxicity collected from various literature sources, USFDA, EMA labels, and other online sources (Bindiganavile et al., 2021;Constable et al., 2022;Davies et al., 1983;Fraunfelder & Fraunfelder, 2021;Gherghel, 2020;Liang et al., 1996;Tehrani et al., 2008), which on data curation and removal of duplicates, reduced to 424 drugs (Ashburn and Thor, 2004;Frederick T. Fraunfelder).For both training and screening datasets, drug features were calculated based on the constitution, topology, charge, and molecular properties using the Biotriangle webserver (http://biotriangle.scbdd.com/chemical/index/).The data pre-processing using the Knime analytics platform (version 4.6.1)removed the features with less than 0.07 variance.Linear correlation was performed for all column pairs (two-tailed) followed by correlation analysis with 0.7 as the correlation threshold which reduced the dimensionality of the features.

Machine learning model development
Machine learning models were developed based on supervised learning algorithms and advanced neural network (ANN) using Python (Ver 3.6).The predictive models were obtained by running the parameters over various supervised machine learning algorithms such as k-Nearest neighbors, Random Forest, a particular class of decision tree C4.5, XG Boost, Support Vector Machines, and Naive Bayes probabilistic techniques.Additionally, k-fold cross-validation (k ¼ 5) was applied to accurately estimate the model's predictive performance.For the implementation of supervised learning models, the Sklearn library was used where substrates were denoted as one and non-substrates as zero in the curated dataset.
The predictions of the test set from six base models (mentioned above) were used as the input data for the logistic regression model to obtain consensus predictions.In addition, ANN was also used to predict the interactions between the drug and hOCT1.It is a particular type of algorithm which consists of 3 layers, namely, the input layer, hidden layer & output layer.In our model implementation, the input layer consisted of 7 neurons and the two hidden layers consisted of 6 neurons, each using a rectifier activation function.
The output layer had one neuron, which used the sigmoid activation function to get the final predictions.
The performance of the developed model was assessed based on several metrics such as accuracy, precision, recall, and root mean square error.The screening dataset (systemic drugs causing ocular toxicity) was run through the developed model, and the molecules were predicted either as substrate or non-substrate for hOCT1.
Since the sequence identity was less than 40% for the target and template protein, the non-template loop regions of homolog were refined with the Prime loop refinement tool using variable dielectric surface generalized born solvation model (VSGB) as solvation model and optimized potential for liquid simulations 4 (OPLS4) force field.Since hOCT1 is a membrane protein, the structural model was oriented according to the orientations of proteins in membranes (OPM) (Lomize et al., 2006).The loop regions of the protein structure containing less than five amino acid residues were refined using the default loop sampling method, six to eleven amino acid residues using the extended loop sampling method, and more than eleven amino acid residues were refined using the ultra-extended loop sampling method.
Further, the protein was prepared using Protein Preparation wizard of Schr€ odinger's Maestro suite to perform molecular docking, MD simulations, and metadynamics.The protein structure was prepared by assigning bond orders, adding missing hydrogen atoms, creating zero-order bonds to metals, creating disulphide bonds, pre-processing, and removing water molecules beyond 5 Å from the Het group.Further hydrogen bonds were assigned using PROPKA at pH 7.4 and energy was minimized for the protein structure.The accuracy of the developed model was evaluated using a protein reliability report and Ramachandran plot, which describes the plot of the torsional angle between phi (u) and psi (w) amino acids in the protein and thereby gives information on the allowed and disallowed conformations of the developed homolog structure (Ramachandran & Sasisekharan, 1968).

Preparation of ligand molecules
Ligand preparation was performed using LigPrep module of Schr€ odinger's Maestro suite.For ligand (known substrates/non-substrates of hOCT1 and systemic drugs causing ocular toxicity) preparation, 3D structures of all the ligand molecules were imported from PubChem.The ligand structures were then prepared at a pH of 7.4 ± 0.2 in the OPLS4 force field of the protein.The prepared ligands were used for molecular docking.

Molecular docking
Molecular docking was performed using Glide module of Schr€ odinger's Maestro suite.Molecular docking was performed to obtain the protein (hOCT1) ligand complex.A receptor grid was generated around the protein structure to define a search space for the docking calculations.A grid of size 30 Å was created around the bound ligand to cover the entire protein space, which was further docked with various prepared ligand molecules using the glide dock function with default parameters unless otherwise mentioned.Extra precision docking (XP) was performed for prepared ligand.The known substrates for OCT1, tetraethyl ammonium (TEA), and 1-methyl-4-phenylpyridinium (MPP) were used as the standards to optimize and validate the model's binding sites.

Molecular dynamic simulations
MD simulations were performed using Desmond module of Schr€ odinger's Maestro suite.We used MD simulations to visualize the interaction and binding stability of the protein and ligand complex (Hollingsworth & Dror, 2018).The protein (hOCT1) and ligand molecules were placed in a system containing TIP3P water molecules, and the membrane region contained 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphocholine (POPC) molecules in a simulation box with buffer size of 10 � 10 � 10 Å.The system was neutralized by adding a calculated number of ions to account for electrostatic neutrality (Gapsys & de Groot, 2020;Hub et al., 2014).The prepared system was subjected to 100 ns MD simulations under the constant number of molecules, pressure, surface tension, and temperature (NPcT).Protein root mean square deviation (RMSD), ligand RMSD, protein root mean square fluctuation (RMSF) and protein-ligand interactions were analyzed at the end of MD simulations.

Metadynamics
Metadynamics was performed using Desmond module of Schr€ odinger's Maestro suite.The translocation of substrates was enhanced by performing metadynamics to visualize the movement of ligands along the protein molecule (hOCT1).MD simulations were performed for 10 ns to equilibrate and relax the system before metadynamics simulation.The distance between the center of mass of the protein molecule and the center of mass of the ligand was used as the collective variable for the metadynamics simulations.A wall-length limit of 35 Å was used to prevent the distancing of ligand molecule from the protein.Since hOCT1 is embedded in the membrane, NPcT ensemble class was used to maintain the constant number of molecules, pressure, surface tension, and temperature.Gaussian height was set at 0.03 kcal/mol with a width of 0.05 Å and applied at intervals of 0.09 ps.The total simulation time for metadynamics was set as 40 ns with a recording interval of trajectory at 40 ps.The temperature was set to 300 K with a pressure of 1 bar and surface tension of 0 Å.The time and wall length were optimized for visualizing the transport process based on the trial and error method.The system's free energy was plotted as the function of the distance moved by the ligand molecule from its initial position during the metadynamics simulation.These energy profile diagrams were used to classify the ligand molecules as substrates or non-substrates of hOCT1.

Dataset preparation
Preparing a dataset with necessary features is the key to a successful model since the physicochemical and structural properties of drugs play a crucial role in their binding with the transporters.Initially, various constitutional, molecular, charge and physicochemical features were obtained for substrate and non-substrates of hOCT1.
An initial 96 features were reduced by variance filter and correlation analysis.Several structural and physicochemical features such as logP, weight, number of sulfur atoms, and others were selected based on their correlation values, a few of which has also been previously reported by us (Baidya et al., 2020).A few features were also selected based on their relevance in substrate translocation, such as topological surface area (TPSA) and the number of acceptor groups (Table S1).Based on the distribution pattern of features among substrate and non-substrate of hOCT1, logP, hydrophilic index, and the number of sulfur atoms are considered significant contributors, followed by TPSA and the number of h-bond donors in the molecules (Figure 1).Most of the molecules found beyond the range of threshold values had a higher probability of being non-substrate.A total of 424 systemic drugs causing ocular toxicity were collected in our database and further categorized based on their therapeutic use (Figure 2).

Machine learning model performance
The prepared dataset has a wide range of physicochemical property values which were normalized before the training.K-fold cross-validation (k ¼ 5) was added to the model, which divides the dataset into k-subsets (or folds) of equal size to test the data on 1-fold and use (k-1) folds to train the model.The process is repeated k-times for validation, each for a different fold.
A supervised machine learning model was developed to predict whether the systemic drugs causing ocular toxicity are a substrate or a non-substrate for hOCT1 based on various features.A total of six supervised learning models were implemented initially to compare their accuracy.Since the training set was small and numeric features were used, there was a possibility of overfitting, which was reduced by using a particular class of Decision trees (C4.5) as one of the base models.All six models were trained with k ¼ 5 as k-fold cross-validation.
Consensus modeling can improve the overall performance of the machine learning model, which comprises several aggregation techniques to obtain better results by combining the strengths of individual algorithms.To overcome the limitations of the individual models, we applied a consensus model based on the concept of stacking and blending, where the predictions from the test dataset of the individual model are the training data for the logistic regression model development.Compared to base models, the weighted consensus approach consistently achieved more favorable values across all evaluation metrics, indicating overall improvements in inaccuracy and stability (Table 1).ANN also classified the drugs as substrate (indicated by 1) and non-substrate (indicated by 0) of hOCT1.The model ran with 100 epochs and presented an accuracy of 80% to 85%.Predictions from the supervised and ANN model classified the drugs as substrate or non-substrate of hOCT1.

Machine learning model predictions and characteristics of predicted OCT1 substrates
A screening dataset of systemic drugs causing ocular toxicity was screened through the developed logistic regression model and ANN model to predict their interactions with OCT1.Out of 424 molecules screened, 125 drug molecules were found to be substrates for OCT1.Since we developed the logistic regression and ANN models, predicted substrates (overlapped) from both models could be the potential substrates of OCT1.Interestingly, these drugs were not reported earlier (n ¼ 125) for any interaction with OCT1 (Table S2) and are found to be the substrate of OCT1 from our analysis.These drugs were further categorized based on their therapeutic use, which shows that more than 43% of anti-infectives (35 out of 80), 28% of CNS (20 out of 71), 21% of CVS (15 out of  Structural, molecular and physicochemical features of substrates (blue) and non-substrates (orange) of organic cation transporter 1 (OCT1).The distribution pattern of structural, physiological and molecular features among substrate and non-substrate of OCT1 demonstrates that log P, hydrophilic index, and the number of sulfur atoms are significant contributors for classification of substrates, followed by topological polar surface area (TPSA) and the number of h-bond donors in the molecules.70), and 12% of anti-neoplastic (8 out of 64) from our dataset were predicted as OCT1 substrates.Computer simulations confirmed these predictions to understand the molecular interactions at the atomic level and visualize the movement of molecules through the transporter.

Development of hOCT1 homolog structure
Computer simulations visualize the interactions between various drugs and hOCT1 at an atomic level and validate the predictions from machine learning models.A homology model was developed for hOCT1 since it has no X-ray crystallographic structure to date (Fiser, 2010).In our studies, GLUT-1 in complex with bound inhibitor (PDB ID-5EQG) was selected as a template based on the identity score, optimum resolution (lesser the value, more the similarity, (2.9 Å)), degree of mutations (none), species similarity (human) and presence of bound ligand.Since the sequence identity between the template (5EQG) and target (OCT1) was less than 40%, energy-based modeling was performed to develop the homolog.The earlier homology model of OCT1 was developed with a template showing only 29% amino acid sequence similarity; however, the quality of the model was evaluated based on the location of substrate-binding amino acids in a single structural epitope and contact with hydrophilic molecules indicating the reliability of the model (Popp et al., 2005).In our studies, these amino acids were located in the core protein region facing the hydrophilic/aqueous medium.Therefore, the developed homology model of OCT1 was consistent and validated using a protein reliability report and Ramachandran plot, as shown in Figure 3.
The prepared model had a few outliers concerning the energy dihedrals and unusual B factors because the model was built with less than 40% sequence identity.However, when the protein reliability report was generated for the region around the ligand within a radius of 10 Å, the outliers were reduced, indicating that the prepared model would not impact the molecular interactions in the binding pocket (Figure 3).The Ramachandran plot showed that the majority (99.27% for hOCT1) of the amino acid residue lies in the favorable and allowed region, indicating the model's suitability for further interaction studies (Ramachandran & Sasisekharan, 1968).The prepared model was energy minimized and used to perform further computational studies.

MD Simulations
MD simulations were performed to equilibrate the docked hOCT1-ligand complex and to visualize the binding stability and interactions of drug-hOCT1.The available literature offers limited predictive value for drug-OCT1 interactions due to structural differences among the substrates (Meyer & Tzvetkov, 2021).Due to the poly-specific nature of OCT1, there is a possibility for the studied drugs to bind at different binding sites and therefore, while performing docking, a grid was generated around the protein rather than a specific site.As reported earlier, the docked structures of TEA and MPP showed interactions with similar amino acids (Phe159, Trp217, Phe244, Asp474) (Koepsell, 2004), validating our homology model and demonstrating that the ligand was docked precisely at the binding pocket of the hOCT1 protein.
The best suitable docked pose of each ligand was selected for MD simulations based on their docking score and interactions with the protein.As these simulations are computationally extensive and require a lot of computational power and time, MD simulations were performed for randomly selected drug molecules from various therapeutic categories.The simulations indicate the ionic and hydrophobic interactions of Trp217 and Asp474 with the quaternary nitrogen of TEA and MPP, as reported earlier (Koepsell, 2004(Koepsell, , 2011;;Meyer & Tzvetkov, 2021).However, in the current study, the template (PDB ID: 5EQG) chosen for the homology model was different than earlier reported, which explains the additional hydrophobic interaction of TEA and MPP with other amino acids, such as Phe244 and Phe159.
MD simulations are known to be independent of the simulation box size when it exceeds a distance of at least 10 Å or three solvation layers from the protein to the box edge, based on which we fixed our simulation box size as 10 Å distance from each side of the protein (Gapsys & de Groot, 2020).Most of the reported literature supports the simulation time of 100 ns to visualize the interactions between ligands and proteins (Amir et al., 2019;Ghosh et al., 2021;Koshy et al., 2010;Schlessinger et al., 2018).Hence 100 ns was chosen as the appropriate time for MD simulations.MD simulations do not simulate molecules' transport direction or entry through the transporter but show molecular interactions between the docked ligand and protein, whether a substrate or non-substrate.Hence, molecular interactions were visualized only between substrates (known and predicted) and hOCT1, as shown in Figure 4. Figure S1 indicates RMSD plots of protein and ligand along with their interactions.The RMSD plots from MD simulations represent the binding stability and equilibration of the ligand-protein complex (drug-hOCT1).A significant fluctuation in the RMSD of protein indicates a conformational change during simulation.In our simulation studies, few deviations were observed with protein RMSD; however, the order changes were within 1 to 3 Å, which is acceptable (Weng et al., 2021).Whereas ligand RMSD as plotted on the secondary Y-axis, indicates the stability of the ligand/drug concerning its binding site (Hermanto et al., 2022).In our studies, the drug molecule's ligand RMSD was stable for most of the simulation time (Figure S1).However, a second RMSD shift was observed after 90 ns for captopril, indicating the drug movement away from the binding pocket (Al-Karmalawy et al., 2021).For other drugs such as cyclophosphamide, risedronate, and sulfadiazine, initial fluctuations in RMSD were observed, however as the simulation progressed, these fluctuations were minimized (Cyclophosphamide �25 ns, Risedronate �10 ns, Sulfadiazine �50 ns).Such fluctuations may arise as the initial docked pose of the ligand might not be the most stable conformation when presented in the solvation medium (Shoichet et al., 1999).MD simulations could enable the alteration of ligand conformation to achieve a stable state, indicated by initial ligand RMSD fluctuations (Liu & Kokubo, 2017).Ligandprotein interactions were used to visualize the interactions between ligand atoms and amino acid residues of hOCT1, such as hydrogen bonds, hydrophobic and ionic interactions, and salt and water bridges (de Freitas & Schapira, 2017).Most of the drug substrates formed hydrophobic interactions with Trp217 (�90%) and Ile446 (�80%), followed by Lys214 (�70%), which formed hydrogen bonds, hydrophobic and ionic interactions.The other highly interacting amino acid residues were found to be Phe159 (�60%), Ser470 (�60%), Tyr221 (�55%), Asp474 (�55%), and Phe244 (�50%).Nevertheless, we were not able to visualize the movement of molecules across the transporter by MD simulations, and therefore the equilibrated structure of drug-hOCT1 was further subjected to metadynamics.Since the conventional MD simulations require a longer duration to simulate the transport movement across the transporter, which is practically not feasible owing to the tremendous amount of computational resources, hardware capacity, and time required for simulating at an atomistic level-metadynamics simulations can be used as an alternative tool to visualize these movements.

Metadynamics simulation
Metadynamics is an atomistic computer simulation tool used to perform accelerated simulations of the biological event by forcing the system across the energy barriers using a series of well-defined Gaussian energy functions to boost the potential energy of the collective variable, thereby speeding up the process (Bussi & Laio, 2020).The metadynamics approach is employed to visualize biological processes like the transport of a ligand molecule through the transporters (ligand-protein equilibrated complex) (Nagy et al., 2021).Apart from improving the simulation timescale, metadynamics also enhances the sampling method by utilizing collective variables whose values directly influence the biological process (Valsson et al., 2016).This study employed metadynamics simulations to classify systemic drugs causing ocular toxicity as substrates and non-substrates of hOCT1.The rationale for selecting drugs (AI predictions) for metadynamics simulation was based on their clinical use and the associated toxicities.To ensure the diversity, the drugs were selected from multiple therapeutic categories such as antiinfectives (6), cardiovascular drugs (9), central nervous system drugs (4), and anti-neoplastics (3).
The known substrates and non-substrates of OCT1 were used to identify the discrete pattern in their associated free energies concerning their movement during the simulations.Upon different trials, it was found that a wall length of 35 Å with a simulation time of 40 ns revealed that around 75% of substrates displayed their lowest energy state close to their initially docked site, and 88% of non-substrates displayed their lowest energy state outside the initially docked site, thereby validating the model with a mean accuracy of 81%.As shown in Figure 5, the energy profile patterns of substrate and non-substrate showed an intriguing pattern as the substrate molecules displayed their lowest energy state very close to their initially docked position, whereas the non-substrates had their lowest energy far away from their initially docked position.This distinct pattern indicates that the substrate molecules were more stable near the binding pocket of the transporter, while non-substrates were stable outside the initial binding pocket.The molecule's position with minimum free energy was used to classify the systemic drugs causing ocular toxicity as substrates and non-substrates of hOCT1 (Video S1).The energy profile graph for various predicted substrates and non-substrates with their classification as substrate or non-substrate is given in Table S3.

Predictions of potential OCT1 substrates from machine learning and computer simulations among systemic drugs causing ocular toxicity
Our earlier studies proved the functional role of OCT1 in the eye in transporting various cationic drugs from systemic circulation to the tear (unpublished data) and across the cornea when administered systemically, intravitreal, and topical (Nirmal J 2010;Nirmal et al., 2013;Velpandian et al., 2012).In contrast, our current study focuses on drug-OCT1 interactions indicating membrane transporters could be a potential portal for the systemic drugs (29% of 424 drugs from our database) into the eye.Some of the OCT1 predictions not known earlier are given in Table 2, along with their therapeutic class and associated ocular toxicities.Several betablockers, including Atenolol, Nadolol, Labetalol, and Pindolol, were found to be substrates for OCT1, as reported earlier using in vitro studies (Guo et al., 2018;Misaka et al., 2016).Interestingly, from our predictions, 33 out of 62 (53%) sulfur-containing drugs were predicted as OCT1 substrates, including busulfan, tamsulosin, dapsone, sulphanilamide, sulfacetamide, and sulphadiazine, indicating the presence of the  sulfur group could be one of the crucial features of OCT1 substrate.Moreover, based on the correlation analysis performed during training dataset preparation, the number of sulfur atoms in the molecular structure was selected as one of the essential features to be identified as OCT1 substrate.However, more detailed invitro and invivo studies are required to confirm this finding.A recent study also reported that sulfur could be an additional factor facilitating the transport of drugs through OCT1 (Redeker et al., 2022).

Conclusion
In our study, we used machine learning, MD simulations, and metadynamics to predict the drug-OCT1 interactions, which could help understand the entry of systemic drugs into the eye-thereby, the ocular toxicity.The findings from our study are: a) predictions from the artificial intelligence model revealed potential OCT1 substrates (n ¼ 125)-not known earlier, b) our predictions demonstrate that the sulfur-moeity in drugs could be an additional factor facilitating the transport of OCT1 substrates, c) metadynamics studies can be used to classify the drugs as substrate or non-substrate based on their free energy concerning their movement during the simulations.Though previous studies have predicted drug-OCT1 interactions using artificial intelligence, to the best of our knowledge, for the first time in the current study, we have also used metadynamics for classifying drugs as substrates and non-substrates for the transporters.This high-throughput screening approach can be further explored to advance the understanding of drug interactions with other transporters.However, these predictions need further validation by invitro and invivo studies to improve our understanding of the entry of systemic drugs into the eye.

Figure 2 .
Figure 2. Therapeutic classification of systemic drugs causing ocular toxicity.The screening dataset consists of 424 systemic drugs causing ocular toxicity which were categorized into different categories based on their therapeutic use.

Figure 1 .
Figure1.Structural, molecular and physicochemical features of substrates (blue) and non-substrates (orange) of organic cation transporter 1 (OCT1).The distribution pattern of structural, physiological and molecular features among substrate and non-substrate of OCT1 demonstrates that log P, hydrophilic index, and the number of sulfur atoms are significant contributors for classification of substrates, followed by topological polar surface area (TPSA) and the number of h-bond donors in the molecules.

Figure 3 .
Figure 3. Evaluation of homology model of organic cation transporter 1 developed using human glucose transporter as template.A and B shows protein reliability report of homology model of protein backbone and region around ligand (10 Å), respectively.Various protein properties are generated based on their structure and conformations, such as steric clashes, bond length and angle deviations, backbone and sidechain dihedrals, planarity, torsions, missing atoms, and stereochemistry with their allowable limits indicating the quality of the developed protein model.C, shows Ramachandran plot of homology model, indicating the location of amino acids in favourable (orange), allowed (yellow), or disallowed (white) regions.It displays the protein dihedrals for all the amino acid residues in the protein where triangle represents glycine, squares represent proline and other amino acid residues are represented by circle.The developed homologue model indicates the presence of glycine (triangle) in the disallowed region, indicating the steric hindrance between the C-beta methylene group of side chains and main chain atoms.However, since glycine has no side chain (absence of methylene group), it does not possess steric hindrance and hence can be accepted in all the quadrants of the plot.

Figure 4 .
Figure 4. Molecular interactions of known substrates with organic cation transporter 1 for a simulation time of 100 ns, evaluated by molecular dynamic simulations.Tetraethyl ammonium forms hydrophobic interactions with Trp 217 and Phe 244, 1-methyl-4-phenylpyrinidium forms hydrophobic interactions with Phe 159, Trp 217 and Phe 244, metformin forms hydrophobic interactions with Phe 159 and Tiotropium forms hydrophobic interactions with Phe 159, and Phe 244.

Figure 5 .
Figure5.Energy profile graph of substrates (TEA, MPP, metformin, Tiotropium) and non-substrates (propranolol, bupropion) of organic cation transporters, obtained through metadynamics simulation studies.TEA: Tetraethyl ammonium, MPP: 1-methyl-4-phenylpyridinium.A distance-based collective variable was used for the metadynamics simulations with a specified center of mass for protein and ligand.A total simulation time for metadynamics was set as 40 ns.The energy profile patterns indicate that the substrate molecules were stable near the binding pocket of the transporter, while non-substrates were stable outside the initial binding pocket.Therefore, the position of molecule where it possesses the minimum free energy was used to classify substrates and non-substrates.

Table 1 .
Evaluation metrics of supervised learning model and artificial neural network.The predictive models were obtained by running the parameters over an artificial neural network (ANN) and various supervised machine learning algorithms such as k-Nearest neighbors (KNN), Random Forest, Decision tree C4.5, Naive Bayes, XG Boost, and Support Vector Machines (SVM).Additionally, k-fold cross-validation (k ¼ 5) was applied to accurately estimate the model's predictive performance.Table 1(a) represents evaluation metrics of different folds as F1, F2, F3, F4, and F5.Table 1(b) represents the overall evaluation metrics of the supervised learning model and ANN.

Table 2 .
Predicted potential organic cation transporters 1 (OCT1).Predictions of OCT1 substrates from machine learning models, further validated by molecular dynamic simulations and metadynamics, revealed drug-OCT1 interactions not known earlier.