{'fraction_non_alphanumeric': 0.03217249796582587, 'fraction_numerical': 0.01708706265256306, 'mean_word_length': 5.529005524861878, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Current Advances in Studying Clinically Relevant Transporters of the Solute Carrier (SLC) Family by Connecting Computational Modeling and Data Science

Organic anion and cation transporting proteins (OATs, OATPs, and OCTs), as well as the Multidrug and Toxin Extrusion (MATE) transporters of the Solute Carrier (SLC) family are playing a pivotal role in the discovery and development of new drugs due to their involvement in drug disposition, drug-drug interactions, adverse drug effects and related toxicity. Computational methods to understand and predict clinically relevant transporter interactions can provide useful guidance at early stages in drug discovery and design, especially if they include contemporary data science approaches. In this review, we summarize the current state-of-the-art of computational approaches for exploring ligand interactions and selectivity for these drug (uptake) transporters. The computational methods discussed here by highlighting interesting examples from the current literature are ranging from semiautomatic data mining and integration, to ligand-based methods (such as quantitative structure-activity relationships, and combinatorial pharmacophore modeling), and finally structure-based methods (such as comparative modeling, molecular docking, and molecular dynamics simulations). We are focusing on promising computational techniques such as fold-recognition methods, proteochemometric modeling or techniques for enhanced sampling of protein conformations used in the context of these ADMET-relevant SLC transporters with a special focus on methods useful for studying ligand selectivity.

SLC family which were identified to have major implications in the pharmacokinetics of drugs. Mostly these transporters are uptake transporters, only MATE1 and MATE2K -the Multidrug and Toxin Extrusion transporters (MATEs; belonging to SLC47A subfamily) -are representatives of efflux pumps that transfer substances out of cells (mainly expressed in kidney and liver cells). The other discussed transporters herein are uptake transporters belonging to the families of Organic Anion Transporting Polypeptides (OATPs; SLCO), Organic Anion Transporters (OATs; SLC22A), Organic Cation Transporters (OCTs; SLC22A), and Organic Carnitine Transporters (OCTNs; SLC22A).
The tissue expression of the different SLC members discussed in this review is quite heterogeneous and was already described extensively by others (summarized in Table 1) [84,97]. However, predominant organs of expression include liver, kidney, and to a lesser extent intestine. Pharmacological barriers in these organs are crucial for absorption, metabolism, distribution, and excretion of drugs and any impairment in the function of one of these transporters can therefore lead to adverse effects or toxicity by altered pharmacokinetics of the drug or food ingredient administered. For instance, treatment with the immunosuppressant cyclosporine can result in statin-induced myopathy when administered at the same time since statins are known to be substrates for OATP1B1. Inhibition of OATP1B1 by cyclosporine therefore leads to increased plasma concentration of different statins [86].
OATP1B1 and OATP1B3 are exclusively expressed in hepatocytes. Due to their high sequence similarity (~80%) they transport common substrates and are inhibited by common inhibitors. Nevertheless, there are some compounds (e.g. pitavastatin for OATP1B1 [41], cholecystokinin octapeptide CCK-8 for OATP1B3 [48]) which are selectively binding to only one of the two transporters with high affinity.
OATP2B1 is ubiquitously expressed, with highest expression levels at basolateral membranes of hepatocytes [97]. Due to its more distant relationship to OATP1B1 and OATP1B3 (~30% sequence identity to both), also substrate and inhibitor profiles are overlapping to a lesser extent [55,113].
Expression of OATP1A2 in liver cells has been observed in cholangiocytes, possibly accounting for the re-uptake of xenobiotics from bile [67]. In addition, common expression of OATP1A2 and OATP2B1 at the luminal membrane of intestinal absorptive cells is potentially implicated in drug-food interactions [28]. Specifically, it has been found that components of fruit juices (especially naringin) cause the inhibition of OATP1A2, which in turn affects oral bioavailability of fexofenadine (substrate of OATP1A2) [7].
OAT1, and OAT3 are expressed at the basolateral membrane of kidney cells, pointing to their concerted interplay in the uptake of substrates and related drug-drug interactions (DDIs) [16,29,44]. For example, it has been found that probenecid inhibits methotrexate transport by OAT1 and OAT3 [88]. OCT2 is expressed at the luminal membrane of kidney cells [37] and its organ of highest expression is liver (rather than kidney).
OCTN1 and OCTN2 are highly expressed in kidney as well [76,124]. OCTN1 and OCTN2 are associated with several pathologies, such as inflammatory bowel disease, primary carnitine deficiency, diabetes, neurological disorders, and cancer [94]. Interestingly, in case of unresectable gastrointestinal stromal tumors treated with imatinib therapy, polymorphisms in OCTN1 and OCTN2 are associated with a prolonged time to progression in GIST patients receiving imatinib therapy [4]. In a different study, OCTN1-mediated uptake of cytarabine into tumor cells in a cohort of acute myeloid leukemia (AML) patients could be related to reduction in the development of resistance to chemotherapy [27].
Human MATE proteins have been shown to be widely distributed across different body tissues, including liver, skeletal muscles, testis, and kidney. Efflux-and uptake-activity of human MATEs have been observed to be pH-dependent [66]. Nevertheless, in renal cells MATE1 serves as an efflux pump, mediating the transport of the substances from kidney into urine. The OCT2-MATE1 interplay in uptaking and effluxing compounds in the kidney shows that also little related transporters can be conjointly involved in clinically relevant DDIs. MATE2-K is the human kidney-specific MATE2 active splice variant of MATE2 [84]. In addition, the ubiquitously expressed variant of MATE, MATE2-B, has also been detected. Nevertheless, this variant probably does not have any functional activity [76].
For being able to assess a compounds' risk to interact with any of these ADMET-relevant transporters, having predictive in silico models for all of them at hand would be useful at an early stage in the drug discovery pipeline. As will be discussed in the following sections, such efforts are complicated by limitations in data availability in the domain of ADMET-relevant SLC transporters. Other challenges include the overlapping substrate or inhibitor specificities across these transporters which are arising from the promiscuous nature of these proteins [111], as well as from the substrate/inhibitor promiscuity.
In this review article, data science and computational modeling approaches for unravelling ligand-transporter interactions for the important class of clinically relevant SLC transporters are discussed. In chapter 2, a general overview on important ligand-and structurebased modeling methods which are traditionally being used in computer-aided drug discovery is provided. Next, challenges arising from data sparseness on both the ligand and protein side in the field of clinically relevant SLC transporters are discussed (chapter 3). Chapter 4 is providing more details about different ligand-based (LB) and structure-based (SB) modeling methods that have been used in the context of clinically relevant SLC transporters and discusses insights that were delivered by the different studies. The emphasis is further put on the combination of LB and SB methods to provide a more comprehensive picture of ligand recognition. In chapter 5, studies and different Kidney, skeletal muscle, brain and placenta Roth et al. [97]. Br J Pharmacol, 165(5), 1260-1287.

OATP2B1
Liver (basolateral membrane), intestine, placenta, heart, mammary gland, brain Roth et al. [97]. Br J Pharmacol, 165(5), 1260-1287. methods providing insights into transporter selectivity are highlighted. Finally, the challenge of including information on transporter flexibility into the SB-modeling procedure is discussed (chapter 6) and some recent developments in the field are highlighted (such as the use of elastic network models for conformational sampling). In conclusion, a computational workflow for studying ligand interactions with clinically relevant SCL transporters is proposed, interconnecting ligand data integration, data curation and analysis, as well as LB and SB modeling techniques in order to come up with binding mode hypotheses.

Molecular Modeling Approaches at a Glance
Computational approaches have become a standard paradigm in area of preclinical drug discovery [78]. Molecular modeling is an interdisciplinary field that incorporates theoretical concepts and efficient computational algorithms to study chemical phenomena [65]. The main idea is to use approximative mathematical models while being able to predict behaviour of chemical systems as closely as possible [129]. In the field of computer-aided drug design (CADD), molecular modeling methods are generally divided into two categories: (1) Ligand-based (LB) and (2) Structure-based (SB) modeling. LB modeling is also known as 'indirect drug design', since in the modeling approach protein information is not taken into account. On the other hand, SB molecular modeling techniques enable to study proteinligand complexes at the molecular level. SB drug design is also termed 'direct drug design' since in these cases 3D information of the target protein or a structural model built on basis of (phylogenetically or structurally) a related protein template is included into the modeling process.
In terms of LB modeling, a system of interest is represented mostly as a statistical model incorporating knowledge on the associated compounds (substrate or inhibitor data) possessing experimentally determined activities against a particular target. Compound series can be used to derive an abstract representation of important molecular features being crucial for a binding event. This approach is known as LB pharmacophore modeling and it is especially useful to overcome difficulties in aligning and searching for structurally little related compounds which persist similar in terms of chemical features (so-called "scaffold-hopping" concept) [45,128]. Another widely used approach is Quantitative Structure-Activity Relationship (QSAR) which is based on correlating physico-chemical properties (or other representations of the chemical structure, such as molecular fingerprints of compounds) to their biological activity [23]. QSAR modeling and related qualitative approaches (such as binary classification models) have currently become especially popular due to the observed increase of the deployment of machine learning and deep learning algorithms in CADD. Especially classification models (models being trained to distinguish compounds into e.g. 'active' and 'inactive' class) can subsequently be used for virtual screening of chemical databases with the aim to detect new (potentially active) compounds. Qualitative models can also aid in the interpretation of ligand profiles which might trigger a compound's activity against a particular target [113]. Beyond the conventional SAR techniques, 3D-QSAR modeling provides a natural extension to the classical QSAR formalism by superimposing 3D ligand structures to retrieve 3D-based ligand features [116].
SB molecular modeling approaches are significantly hampered by the limited number of experimentally resolved structures for membrane transporters. So far (updated on January 2019), the number of the available crystal structures for membrane proteins in Protein Data Bank (PDB) did not exceed~3,5% of all deposited entries. Such a small fraction of resolved membrane proteins is primarily caused by problems when overexpressing membrane proteins in bacteria [101], as well as by the obstacles accompanying crystallization of membrane proteins [89]. These include, for example, finding optimal conditions for crystallization [87], as well as difficulties to account for a membrane-like environment which in turn is being crucial for crystallizing the native form of membrane proteins. However, substantial progress has been made by e.g. improving protein engineering to increase the stability of a protein of interest [110]. The latter issue is commonly solved by using micelle detergents to solubilize membrane proteins for purification purposes [106]. Except for X-ray crystallization, solid-state NMR techniques can be used to resolve structure of membrane proteins [63]. In this methodology, membrane nanodiscs have successfully been applied as membrane-mimetics [98]. If the crystal, NMR, or recently also Cryo-EM structure of a respective target is not available, 3D structural models can be built "from scratch" upon basic thermodynamic principles [39]. These methods are known as ab initio (de novo) structural predictions and are rather rarely used, mainly due to the challenges to efficiently sample the whole conformational space at a feasible time scale [51]. In most of the cases, 3D structural models are created on basis of sequence-homologous proteins with known structure ("homology modeling" [52]), or more evolutionary distant proteins with conserved fold ("fold-recognition methods" or "protein threading" [75]). Fitting of a target sequence onto the 3D coordinates of a sequentially-or structurally-related protein ("template") is further accompanied by a global geometry optimization to satisfy structural restraints imposed by internal coordinates, as well as local optimization to remove steric clashes and/or reduce the noise caused by poorly modeled side chain rotamers [62,99,119]. This procedure is commonly applicable for the transmembrane core of membrane proteins. However, it might not be sufficient for modeling intra-and extra-cellular domains, which are mostly consisting of large loop regions [32]. Loop modeling is a non-trivial step in protein structure prediction, since the loops are intrinsically disordered regions, often requiring enhanced conformational sampling [33]. In addition, loops are usually indeterminate parts in the crystallization process due to their high conformational flexibility and thus low electron density in X-ray diffraction patterns. Homology modeling is subsequently limited by incomplete sequence alignment because of missing loop regions. Modeling of extremely short segments (b3 amino acids) can be satisfied by geometrical constraints of their bond lengths and angles. Template-based loop construction by using a database of known structural fragments (not necessarily with identical sequence) is a common modeling approach for medium-size loops (~10 amino acids). In general, energy minimization combined with MD and simulated annealing is advisable to refine modeled loops. For longer loops (~25 amino acids), de novo coarse grained modeling methods (employing e.g. united residue models) have successfully been applied [49].
Furthermore, molecular dynamics (MD) simulations with enhanced sampling techniques can additionally be integrated to the 3D model building procedure especially to refine low-confidence regions, such as flexible loops (as discussed above) and less structurally-conserved parts of the protein [34]. When performing MD simulations on membrane proteins, one should also account for the substantial role of phospholipid membranes which spatially restrain a protein's 3D structure [82]. The computationally less demanding approach is to treat membrane environment as a mean-field continuum model which replaces explicit protein-lipid interactions by effective interactions being a priori included in force-field parameters of a membrane-protein system [31]. These approaches have become particularly useful for the simulations of large time scale events, such as protein folding, albeit for the cost of lacking all-atom representation of a simulated system. The lack of high-resolution accuracy by using implicit membrane models can be corrected by the explicit representation of lipid bilayers in the simulation box. However, running MD simulations with explicit lipid bilayers might be unfeasible for biologically relevant time scales.
Quality assessment of 3D protein models is required to distinguish the native protein structure from the physically non-relevant states [40]. For this purpose, several scoring functions, including statisticsbased, knowledge-based, physics-based, or their combinations, have been developed. To give an example, ProQM is a statistics-(learning-) based method using Support Vector Machines (SVM) models which are trained on the known structures to predict correct structural features of membrane proteins, such as membrane topology or conserved structural motifs [95]. Another example is C-score which estimates the confidence of target-template alignments based on the foldrecognition methods [127]. TM-score function is a metric of 3D similarity between two proteins when performing structural alignment [131]. The next step in SB modeling is to apply docking algorithms to iteratively search for preferred orientations of a ligand molecule relatively to the protein binding site(s) [68]. Subsequently, protein-ligand complexes can be ranked on basis of scoring functions (knowledge-based [85], physics/force field-based [35], machine-learning [2], and/or empirical [12]), estimating the likelihood of all possible binding poses with respect to the energetically (un)favorable intermolecular interactions. Molecular docking can be performed by using different strategies, ranging from rigid docking, where the protein structure is kept fixed and only the ligand's conformational space is explored [9], to induced-fit docking, where the protein local backbone movements are allowed to adjust the proper ligand-protein binding pose [117]. More sophisticated and computationally demanding methods are treating the whole protein structure as flexible [123]. Docking screens are typically evaluated on basis of score accuracy (by comparing the predicted binding affinities to the experimental ones, if available [8]), enrichment factor (by checking if the docking screens are able to discriminate between known actives and known inactives/decoys [46]), or on basis of prospective validation (by measuring e.g. IC50 values [47]). Traditional docking approaches can be complemented by more accurate free energy calculations not only to estimate absolute free energy of binding [24], but also to study ligand selectivity profiles [3]. Approximations of binding free energies are given by methods like MM/GBSA (Molecular Mechanics/Generalized Born Surface Area), as originally described by Kollman et al. [59]. However, free energy perturbation (FEP) provides better estimates of free energy of ligand binding, and it offers the possibility to directly evaluate the impact of mutated functional groups of the ligand from the energetic point of view [118]. On the protein side, biochemical mutational studies can be informed by docking exercises and vice versa as demonstrated for the interleukin 8-gene [22].
Following the same strategy as in case of LB pharmacophores, SB pharmacophore models can be created upon projection of the important pharmacophoric features of a target-ligand complex to an abstract representation [107]. In the recent past, new techniques combining LB and SB approaches have become popular. For example, proteochemometric modeling (PCM) can outperform traditional QSAR modeling by simultaneous evaluation of the similarity of ligands and targets [115]. The great advantage here is the possibility to extrapolate the activities of known ligands against known targets to novel targets without knowing their 3D structure.

Data Sparseness as a Major Challenge in Computer-aided Drug Discovery & Implications to the Exploration of Uptake Transporters
Chapter2 provides on overview of main computational techniques which are traditionally being used in the drug discovery pipeline. When it comes to investigations on clinically relevant SLC transporters, however, the direct application of above-mentioned methods is far from being trivial. The main obstacle hampering computational studies is caused by data sparseness on both ligand and protein levels in the domain of uptake transporters. In addition, the promiscuous nature of uptake transporters being able to bind both endogenous compounds and xenobiotics [56], considerable transporter flexibility which accompanies translocation processes (such as "rocker-switch" mechanism as proposed for Major Facilitator Superfamily members) [120], as well as the likely existence of multiple binding sites [43,70,74,125], turns all modeling efforts into even more challenging tasks.
From a ligand's perspective, modeling studies are strongly impeded by the inconsistent and mostly insufficient number of high-quality bioactivity data which is spread over different data sources in the open domain. For LB studies on uptake transporters, single-point percentage inhibition data has often become the only source for QSAR modeling [1,14,19,55,58,60,122]. Such models trained on single-point inhibition data can achieve high accuracies, particularly if the measurements were retrieved following the same experimental protocols, as demonstrated by e.g. the classification models by Karlgren et al. [55] (with accuracies between 73% and 92% for the different models). However, even in the case of almost identical assay protocols, other parameters such as the substrate concentration might to a significant extent change the final value of the read out (e.g. percentage inhibition value). In case of the studies by Ahlin et al. [1] and Chen et al. [19] which both measured reuptake inhibition for OCT1 inhibitors by using the same probe substrate (4-4-dimethylaminostyryl-N-methylpyridinium, abbreviated as ASP + ), a different substrate concentration (2 μM in Chen et al. and 1 μM in Ahlin et al., respectively) as well as inhibitor concentration (20 μM in Chen et al. and 50 μM in Ahlin et al.) led to a percentage of around 14% of all overlapping compounds (measured in both papers) being differently classified by the different studies (by using a cutoff of 50%). Obviously, the study by Chen et al. was much more rigorous in assigning the label "inhibitor", since many conflicting annotations with Ahlin et al. turned out to be rather false negatives (data not shown).
In most of these cases, the activity cut-off for binary classification into actives and inactives was set to 50% [1,19,55,60,122]. In some studies, however, the authors applied more stringent activity criteria for setting a treshold by removing weak actives from the data set, e.g. in the study by van de Steeg et al. [114]. Here, the activity cut-off for OATP1B1 inhibitors was defined as ≥ 60%, while OATP1B1 noninhibitors were defined as ≤ 40%. It concludes that all "grey zone" data points (the weak actives) in the range (40;60) were excluded from the dataset. On one hand, removal of weak actives can be beneficial to reduce the noise caused by the variations in the experimental measurements [114]. On the other hand, complete exclusion of weak actives might decrease the applicability domain of such models, e.g. when attempting to predict the activity of compounds which are structurally closely related to those from the "grey zone". For establishing LB pharmacophore models for MATE1 [126] and OCT2 inhibitors [125], classification of inhibition/substrate data into binary classes was done via manual literature searches with the aim to extract recommended activity thresholds proposed by the authors of primary literature sources. For deciding upon the cutoffs for percentage inhibition data in order to generate predictive classification models for OATP1B1, OATP1B3, and OATP2B1 inhibition from diverse data source we have chosen the same strategy [113].
Structurally and pharmacologically distinct compounds can also be partitioned into discrete clusters based on their molecular features, as applied for e.g. the classification of OCT2 inhibitors/non-inhibitors [58]. After the clustering step, different activity cut-offs (inhibitory effect ≥ 50% or ≥75%) were probed to examine which clusters contained the highest fraction of the inhibitors based on the applied cut-off. Specifically, by applying the more stringent cut-off (≥ 75%) the inhibition patterns of identified clusters became more pronounced. Distinct OCT2 clusters were subsequently used for deriving several independent SAR analyses to explain complementary inhibitory mechanisms of OCT2 inhibitors.
Full dose-response curve data has been exploited to much lesser extent. Examples include IC50 measurements for LB studies on MATE1, MATE2-K [5], and OCT2 [122], Km values for OATP1B1 substrate pharmacophore models [17], as well as Ki values for pharmacophore modeling of OCT2 stereoselective binding [79].
If the amount of bioactivity data is not sufficient for model development, direct usage of categorical annotations, such as "substrate", "non-substrate", "inhibitor", or "non-inhibitor", can be applied. For this purpose, Drugbank (containing a collection of marketed or FDA-approved drugs [121]), or Metrabase (primarily containing substrates for OCT1, OATP1A2, OATP1B1, OATP1B3, and OATP2B1 [72]) can serve as rich, open sources for LB modeling. For example, substrate annotations from Metrabase were retrieved for QSAR and PCM modeling to predict OCT1, OATP1A2, OATP1B1, OATP1B3, and OATP2B1 substrates [104]. It has to be pointed out, that using manual activity annotations for binary classification modeling seems to be quite error-prone since it is not clear how data curators decided upon assignment of activity annotations in certain cases. As demonstrated in our current LB studies on hepatic OATPs [113], an extensive comparison of Metrabase annotations for OATP1B1, OATP1B3, and OATP2B1 (non-)substrates and (non-)inhibitors with numerical bioactivity measurements from CHEMBL revealed activity misclassifications for Metrabase data up to 74%. On the other hand, categorical annotations can still be utilized in developing accurate computational models when e.g. performing selective fusion of more independent classifiers and therefore increasing the confidence in the model's predictability [104].
With increasing efforts of making bioactivity data publicly available to the scientific community, new challenges are arising. Nowadays, it is no longer only access to data but the proper usage of data which can provide a competitive advantage in drug discovery. Finally, filtering out high-quality data is essential as well as making use of the possibilities to interconnect different types of data, including data originating from diverse sources. In the light of those efforts, data analytics platforms for creating automated data workflows, such as the Konstanz Information Miner (KNIME [11]) or Pipeline Pilot [108], have become particularly useful. Notably, integrative data mining (i.e., data fusion from different sources) can enrich existing data sets by not only the size in enumerated compounds, but also by obtaining novel scaffolds which can lead to an expansion of the available chemical space as demonstrated recently for hepatic OATPs [113]. Moreover, merging data from multiple independent bioactivity measurements (Km, IC50, EC50, Ki, percentage inhibition data) for a single compound can significantly increase the confidence of bioactivity data. When it comes to QSAR modeling, it is usually not recommended to mix compound data originating from different bioactivity end-points [53,61]. On the other hand, binary classification (e.g. separating substrates from non-substrates) should be independent from a certain assay or experimental protocol used [80]. Another benefit when considering multiple bioactivity measurements is the ability to rationally decide upon activity thresholds for binary label assignment (e.g. active/inactive) by studying the distribution of bioactivities within the data set [113].
As recently demonstrated [113], above mentioned pipelining tools are quite handy for semi-automatically fetching ligand data from different open data sources, such as CHEMBL [10], UCSF-FDA Transportal [83], IUPHAR [90], and Drugbank [121] and Metrabase [72]. It has to be emphasized that ligand data originating from different sources might be highly inconsistent with respect to their structural format used. To overcome this issue, applying a standardization protocol, as e.g. by Atkinson (available at https://wwwdev.ebi.ac.uk/chembl/ extra/francis/standardiser/), has proven to be useful for mapping data from different sources.
From Table 2 and Table 3 it becomes clear that the most comprehensive data set is currently available for OATP1B1 and OATP1B3, with CHEMBL being detected as the most prominent source (1993 and 1972 unique compounds, respectively). This trend is reflected by the high number of LB modeling studies with a special focus on OATP1B1: QSAR/classification models [6,54,55,60,104,105,113,114], PCM [14,104], and LB pharmacophore model [17].
In addition to this quite restricted availability of compound bioactivity data for uptake transporters which limits over all the applicability domain of respective LB models, SB modeling efforts are strongly impeded by the lack of crystal structures in this domain. Only MATE1 bacterial templates are available (NorM from Vibrio cholera21, pdb id: 3mkt; pfMATE from Pyrococcus furiosus, pdb id: 3vvn). Both of the available crystal structures share approximately 24% identity with human MATE1. Except for this case, there are no available homologs for the other clinically relevant SLC transporters. These drawbacks can be diminished by the use of fold-recognition methods to search for structurally related analogues with conserved fold. The rationale behind is that a protein's secondary structure should have been conserved to a higher extent during evolution than its amino acid sequence [52]. An overview of available crystal structures which were predicted as potentially useful templates for SB modeling of ADMET relevant SLC transporters is provided in Fig. 1. Predictions were performed by using pGenThreader prediction server [71]. It appears interesting that all predicted templates (except for chain L from respiratory complex I) belong  to the Major Facilitator Superfamily (MFS) [92]. All protein structures are built up of twelve transmembrane helices. It is interesting to note that sequence similarity among discussed transporters and the detected template structures is generally rather low (10-25%), which obviously slowed down the generation of comparative models on basis of these templates compared to other protein families (numbers of published studies per year including comparative modeling for clinically relevant SLC transporters is depicted in Fig. 2). As visible from Fig. 1

Linking Ligand-and Structure-Based Modeling Methods for Studying Uptake Transporters
Natural increase in available ligand bioactivity data and structural templates also led to an increase of published computational studies on clinically relevant SLC transporters over the years. Fig. 2 shows the time evolution of different computational methods that have been used in the context of uptake transporter modeling. In general, we can observe an increase in in silico studies over the years for both LB and SB approaches. Both approaches have been used early on and we can observe peaks in the emergence of QSAR/classification models and new homology models in 2011/12. In 2017, one publication reported the establishment of new comparative models for OCT1, OCT2, OCT3, OCTN1, and OCTN2, which is also visible as a peak from Fig. 2. Interestingly, MD simulations were integrated only twice into the process of clinically relevant SLC transporter modeling, namely for OAT1 and MATE1 (in 2011/12). With the recent increase in available comparative models for uptake transporters, a revival of MD-based methods can be expected especially for performing enhanced conformational sampling of distinct transporter states in the future (see also chapter 5). Naturally, studies including PCM modeling appeared later in literature, since they incorporate information from the ligand and protein side. Thus, for such a technique to be effective it requires a minimum amount of ligand data and ideally knowledge on potential protein binding sites.
Also, from a static perspective, LB and SB methods have been used to approximately the same extend (e.g. 26 studies on QSAR/classification modeling; 26 studies including comparative modeling). If a fair amount of compound bioactivity data for a particular transporter is available (at least a few hundred unique compounds), LB approaches can deliver quicker and more comprehensive insights into important molecular features driving compound affinity towards a particular target than SB approaches. On the other hand, traditional SAR-based methods do not account for the hypothetical presence of distinct binding sites. Such drawbacks can be diminished by sub-structural pattern analyses such as studying pharmacological profiles of enriched scaffolds in the data set. The hereby retrieved congeneric SAR series are likely to interact with the same binding pocket and differences in pharmacological activities are likely triggered by subtle modifications at the side chain level [113]. Thus, such compound series are useful collections for subsequent structure-based docking studies.
In addition, LB pharmacophore models can be useful in detecting important pharmacophoric ligand features which can complement SB docking studies into comparative models. An interesting pharmacophore-based modeling approach where in vitro and in silico ('IVIS') methods are combined has been adopted by Diao et al. to map pharmacophoric features of human OCTN2 inhibitors [25,26]. This methodology has subsequently been used also by Astorga et al. to study the inhibitory profiles of human MATE1 [5]. Specifically, the IVIS procedure aims to iteratively build 3D-pharmacophore models, further used for database screening and subsequent experimental testing of new hits. Afterwards, high-affinity detected compounds are used to rebuild initial pharmacophore hypotheses to perform another round of database screening, and so on. For human MATE1 inhibitors, a common-features pharmacophore has initially been developed by merging pharmacophoric features of both high-and low-affinity MATE1 inhibitors [25]. The aim of mixing high-and low-affinity compounds for building a single pharmacophore model is to detect the minimal essential properties which are crucial for the effective interaction with MATE1. The iterative procedure has finally led to a pharmacophore with two hydrophobic features, one H-bond acceptor and one ionizable (cationic) feature. It is noteworthy to mention that the generation of quantitative pharmacophores has been strongly dependent on the probe substrates used in the in vitro measurements. Specifically, when using 4-4dimethylaminostyryl-N-methyl-pyridinium instead of 1-methyl-4phenylpyridinium, the final pharmacophore model was composed of three hydrophobic features, two H-bond acceptors, and three excluded volumes, spatially arranged in strikingly different configuration from the original one. These findings suggest that human MATEs might contain multiple substrate (and potentially inhibitor) binding sites. Therefore, for developing a complete computational model of MATE inhibitory mechanisms it would be required to build up multiple pharmacophore models per distinct binding site. For OCTN2 inhibitors, the common pharmacophore model revealed three hydrophobic  features and one positive ionizable feature [26]. This IVIS-based pharmacophore hypothesis was partially confirmed by the upcoming study on OCTN2, showing that two hydrophobic features, one Hbond donor, as well as one positively ionizable feature are likely driving OCTN2 inhibitory activity [25]. As already introduced in chapter 3, SB modeling is complicated by the lack of native family member templates. Exclusively human and rabbit MATE1 structural models have been constructed on basis of sequence-homologous templates retrieved by the BLAST algorithm [126,130]. The first structural model for human MATE1 has been built upon the sequence similarity with NorM crystal structure from Vibrio cholerae (pdb id: 3mkt, 35.6% sequence similarity) [130]. After structural model generation, molecular dynamics (MD) simulation has been performed to test the stability of generated human/rabbit MATE1 structural model. In parallel, MD simulations have been performed for the NorM template structure to see if the conformational dynamics of Norm crystal and derived homology models remains conserved. A 50 ns long production MD confirmed overall stability. Interestingly, several helices in the human MATE1 homology model (for example, TM6 and TM9) reoriented and assumed opposite tilt angles when compared to the template. Obviously, a partial closing of the translocation pore was happening for the homology model. This interesting use case highlights the fact that each transporter ortholog might possess its internal dynamics, which cannot always easily be captured by the template structure.
Other reported SB studies rely on the structural similarity with MFS members. For details describing studies for a respective uptake transporter, see Table 4. In the recent past, a multiscale approach for 3D structural modeling of seven human OCTs (OCT1, OCT2, OCT3, OCTN1, OCTN2, OCT6, and FLIPT1) has been published. This study [21]. This study [144] introduces a comprehensive modeling pipeline for tertiary structure prediction, starting from the comparative sequence alignment, and fold-recognition 3D model building combined with ab-initio modeling performed via I-TASSER [145]). In addition, post-translation modifications of functionally relevant structural motifs (e.g. phosphorylation, ubiquitination, and/or glycosylation sites) were predicted via bioinformatic tools (PhosphoSitePlus available at http://www.phosphosite.org/ and NetNGlyc 1.0 server available at http://www.cbs.dtu.dk/services/NetNGlyc/). An integrative approach for 3D structure prediction is followed by the comprehensive evaluation of the modeled structures based on different metrics, including sequence identity between the target and template, query coverage, and consensus Z-score of the top threading programs. Structural analysis of generated models has revealed that the 3D structural models generated here share structural similarity with the human glucose transporter GLUT3 (pdb id: 5c65). Visual inspection of the obtained models further implies 2-pseudofold symmetry, as well as the hypothesis about two distinct functional states (inward-and outward-open).
These observations were fully supported by the structural superposition of 3D generated models of OCTs with the GLUT3 transporter. This use case again demonstrates that phylogenetically unrelated transporters (since OCTs belong to SLC22A and GLUTs belong to SLC2A subfamily) can share the same fold.
Two very promising ways of integrating ligand and protein information are so-called combinatorial pharmacophores as well as PCM modeling approaches. Combinatorial pharmacophore modeling was first reported for OCT2 inhibitors in 2013 [125]. In general, this approach represents a multi-step combinatorial scheme to generate a set of diverse LB pharmacophore models including the available information about their binding modes: First, generated pharmacophore models with identical pharmacophoric features in a close spatial arrangement are grouped in order to reduce the large pool of potential hypotheses and a combinatorial approach is employed to test all possible combinations of different pharmacophore hypotheses. The main idea behind the combinatorial pharmacophores is to study how different pharmacophoric patterns are corresponding to (potential) multiple binding mode hypotheses of uptake transporters. For this purpose, different sub-categories of reference inhibitors are being used -(1) competitive inhibitors (i.e., binding to orthosteric binding site), (2) occluding inhibitors (i.e., noncompetitive inhibitors, occluding the substrate binding site and locking the conformation transformation of the tarnsporter), and (3) allosteric inhibitors (i.e., modulating the transporter's function by binding to the differentallostericbinding site). A use case on MATE1-OCT2 selectivity profiling is presented in chapter 5.
Proteochemometric modeling (PCM) is conceptualized as an advanced extension to the conventional QSAR-based modeling by simultaneous considerations of the similarity between multiple ligands and multiple targets [115]. PCM modeling can thus be categorized as a method at the interface between ligand-and structure-based modeling. The two-dimensional structural sequence information can be integrated into the PCM model either as a whole amino acid sequence, or the preselection of key residues (e.g. those occurring in the binding pocket and/ or other conserved residues) can be performed. ( [64,93]; [138]). Protein sequences can then be reduced to a more abstract representation by calculating the Z-scales which are corresponding to the principal components of multi-property matrices combining different physicochemical properties, such as lipophilicity, volume, and polarity for respective residues [64,115]. In case of uptake transporters discussed herein, PCM modeling was used for the investigations of structural determinants between OATP1B1 and OATP1B3 [14], as well as in a recent study by Shaikh et al. for investigating transporter substrates of OCT1, OATP1A2, OATP1B1, OATP1B3, and OATP2B1 [104].

Selectivity Profiling: Linking Knowledge of Related Uptake Transporters
Since uptake transporters are often co-expressed at pharmacological barriers and generally transport a wide variety of pharmaceutical agents, it is of medical interest to increase the understanding of the interplay of such related transporters. A prominent example are hepatic OATPs -OATP1B1, OATP1B3, and OATP2B1which are responsible for e.g. bile acids uptake (such as taurocholic acid), but also pharmaceuticals, hormones etc. into hepatocytes ( [139], [140]). It is only insufficiently understood to date, how ligand activity and selectivity towards one of the three transporters is achieved. Such knowledge could not only pave the way for functional studies on these transporters (by the use of truly selective tool compounds), but would also increase our knowledge on critical compound/drug properties associated with the onset of clinically relevant drug-drug interactions.
On principle, studies on determining factors for selectivity can include knowledge from the ligand side (QSAR/classification modeling, pharmacophore modeling), the protein side (comparative modeling, molecular docking, virtual screening, MD simulations), or both (PCM modeling or combinations of the latter approaches).
In case of OATP1B1, OATP1B3, and OATP2B1, comparative QSAR modeling has been performed by Karlgren et al. already in 2012 [55] which identified important molecular features for general OATP inhibition (vs. non-inhibition): higher lipophilicity, molecular weight and polarity. Just recently our group identified additional features discriminating hepatic OATP inhibitors from non-inhibitors: higher polarizability, molecular refractivity (corresponding to the distribution of charge over a molecule's surface), and flexibility (expressed as a higher number of rotatable bonds) [146]. Development of in silico models for individual OATP transporters by Karlgren et al. [55] has identified certain differences between the OATP1B and OATP2B subfamily. In contrast to OATP1B1 and OATP1B3, OATP2B1 inhibitory activity has been negatively correlated with nonpolar-and total-surface area, proposing that OATP2B1 inhibitors might be less dependent on polarity than OATP1B1 and OATP1B3 inhibitors. In addition, in our current study we could highlight additional properties to be responsible for OATP1B1 and OATP1B3 versus OATP2B1 inhibition (the latter seem to be more planar, whereas OATP1B members tend to possess a large number of amide bonds) [113].
Another way to explore ligand (and potentially selectivity) profiles is an enrichment analysis in substructures among actives of one target of interest vs the other(s). Again, for hepatic OATPs, this methodology led to a list of enriched scaffolds possessing a certain activity profile (i.e.,OATP1B1 selective inhibition, OATP1B1/OATP1B3 dual inhibition, OATP1B1/OATP1B3/OATP2B1 pan-inhibition). As an outcome, e.g. the pravastatin-like scaffold showed a preferential inhibitory activity for OATP1B1 (over OATP1B3 and OATP2B1) and the cyclosporine-like scaffold accounted for OATP1B1/OATP1B3 dual inhibition. Interestingly, the steroidal scaffold has been found to be enriched in the actives of all three hepatic OATPs. Here depending on the side-chain variations, preferred activity towards one of the targets might be achieved [113].
Pharmacophore modeling is also interesting for studying ligand selectivity across different species. To give an example, human-and rabbit-OCT2 pharmacophore models indicate that despite the similarity of most of the pharmacophoric features (reflected by 83% sequence identity of these two OCT2 variants), there is a difference in the spatial arrangement of hydrogen bonding features [109].
Finally, even ligand profiles and selectivity among uptake transporters of different families might be of interest, in particular if they are commonly expressed at the same pharmacological barrier. A way to tackle this is comparing pharmacophore hypotheses generated for the two targets of interest, like in the case of OCT2 and MATE1, which both are playing a significant role in renal disposition and toxicity (König et al. [141). It has been shown that charge distribution was one of the important factors, favoring the inhibitory activity of one transporter with respect to the other. Specifically, OCT2 inhibitors comprise both zwitterionic and basic functional groups, whereas MATE1 inhibitors are less enriched with basic groups and do not necessarily contain zwitterionic groups [122].
An even more comprehensive understanding of OCT2-MATE1 selectivity profiling was delivered by a combinatorial pharmacophore-based approach [125,126], as introduced in chapter 4. Since MATE1 and OCT2 are commonly expressed in the kidney, it is interesting to learn about their interplay and selectivity switches to better understand transporter-mediated drug distribution and drug elimination processes (König et al., 2011). A combinatorial pharmacophore model approach developed for both OCT2 [125] and MATE1 [126] can therefore reveal which features are shared and which ones are unique for just one of these two transporters. The latter can give hints for selectivity switches at the ligand level. For OCT2, combinatorial pharmacophore modeling has revealed four distinct pharmacophoric hypotheses [125]. An aromatic feature was included in all four hypotheses, suggesting the essential role of pi-pi interactions in the OCT2-ligand recognition. In addition, a cationic charge has appeared in three out of four pharmacophoric hypotheses, which corresponds to previous findings [122]. Xu et al. also compared the molecular weights for the inhibitors matching different pharmacophoric features which provided additional insights into the constitution and/or size of distinctive binding site(s) within the transporter. Following the same strategy as for OCT2, the authors studied multiple inhibitory mechanisms of MATE1 ligands in a follow-up paper [126]. The model reveals significant importance of aromatic rings, as well as hydrophobicity to induce MATE1 inhibition. When compared to the combinatorial model for OCT2 inhibition, it becomes obvious that one of the hypotheses was the same in both transporters, thus proposing one common binding mode hypothesis which can accommodate a substantial number of dual MATE1 and OCT2 inhibitors.
In the future more sophisticated methods, such as multi-label classification might come into play, depending on the availability of compound data with consistent bioactivity measurements for targets under study. Such methods were recently used for studying selectivity profiles of ABC transporters [81].
In addition to the above discussed LB approaches for studying selectivity, the molecular basis for selectivity is delivered by the protein structure, and more specifically by subtle differences in residues interacting with the ligand during binding and transport. To give an illustrative example, attempts to understand selectivity among hepatic OATPs at a molecular level are discussed. Transmembrane regions for both OATP1B3 and OATP2B1 were built on basis of templates from the MFS family, namely glycerol-3-phosphate (pdb id: 1pw4) and lactose permease (pdb id: 1pv6) from Escherichia coli [77]. The model quality has subsequently been validated via docking of the cardiac glycoside digoxin into the putative translocation channel. Based on the 3D models and multiple sequence alignment across OATP family members, the analyses suggest that the pore-facing residue R181 might contribute to the substrate specificity of OATP1B transporters, as this residue is fully conserved across the OATP1B family. In analogy, H579 is hypothesized to be crucial for binding of ligands to members of the OATP2B family and it is found at a spatially adjacent position to R181 [77]. Another SB modeling effort for understanding commonalities and differences of the more closely related hepatic transporters OATP1B1 and OATP1B3 (~80% sequence identity) led to the construction of a series of chimeric proteins between OATP1B3 and 1B1 [38]. The aim here was the determination of structural domains and/or residues responsible for substrate selectivity of OATP1B3, specifically for CCK-8. Homology modeling and molecular docking led to binding mode hypotheses which were further validated experimentally. When replacing TM10 in OATP1B3 with TM10 of OATP1B1 a dramatically reduced degree of CCK-8 transport was observed, indicating that TM10 is indeed crucial for recognition and/or translocation of CCK-8. Using site-directed mutagenesis, key residues for substrate binding namely, Y537, S545, and T550 in TM10 were identified [38].
Using ligand and protein information in conjunction for selectivity profiling can be conducted by using PCM modeling. This technique outperforms conventional QSAR models and can be used to virtually screen for selective compounds that are solely active on a single member of a subfamily of targets [115]. In the field of clinically relevant SLC transporters, PCM modeling was first undertaken for investigating chemical features favoring OATP1B1 and OATP1B3 inhibition [14]. Performing multiple sequence alignment for OATP1B1, OATP1B3, OATP2B1, and OATP1A2 aided in identifying most conserved regions. Further, critical protein residues were prioritized on basis of SB modeling studies previously done for hepatic OATPs [73,77]. This step again demonstrates the usefulness of combining different computational approaches. PCM models were developed for OATP1B inhibition (2-class classification model, i.e., 'OATP1B inhibitor' or 'OATP1B non-inhibitor'), and individual OATP1B1/1B3 inhibition (4-class classification model, i.e., 'OATP1B dual inhibitor', 'OATP1B dual non-inhibitor', 'OATP1B1 selective inhibitor', or 'OATP1B3 selective inhibitor'). When looking at protein properties, only limited conclusions could be deduced from this study. The limited interpretability of target information is caused by the fact that only two proteins were included into the PCM modeling (van Westen et al., 2012). As a future perspective, the authors suggest to apply the PCM procedure to more OATP proteins, as there are bioactivity measurements for 22 OATP isoforms available in CHEMBL. In conclusion, the developed 2-and 4-class classification models united the important molecular features reported from previous computational studies [6,17,55,105], but also provided new information about OATP1B1/1B3 inhibition (e.g. high number of ring bonds).
Pharmacological profiles can also be investigated for a whole group of simultaneously expressed transporters located at the same pharmacological barrier. For example, Shaikh et al. combined QSAR and PCM modeling to perform an extensive exploration of substrate interactions for 13 clinically relevant efflux and uptake transporters, including OCT1, OATP1B1, OATP1B3, OATP2B1, and OATP1A2 [104]. The motivation behind such an extensive study is to predict transport across major pharmacological barriers as a whole by a sequence of computational models. Thus, consensus models by applying various machine learning techniques were finally constructed and the developed models were used to predict the substrate propensity of compounds for blood-brain barrier (BBB) transport.

Accounting for Transporter Flexibility
Since substrate translocation via transporters requires the protein to adapt different conformational states, protein flexibility has to be taken into account also when performing SB modeling [30]. So far, most of the docking studies on uptake transporters were performed by considering only a single transporter conformation (e.g., OAT1 [93], OCT1 [13], OATP1B1 [69], OATP1B3 [38]). However, building a structural model on basis of a single conformation can inherently bias subsequent docking screens since e.g. some ligands might not fit into a narrow binding pocket of a particular conformation. Khuri [57]. In general, multiple comparative models can together provide a more comprehensive representation of a ligand binding event [100]. However, combining drastically different conformational states (such as outward-open with inward) is risky, since allosteric modulations rather happen upon subtle changes in protein conformation and slight rotameric variations of side chains [96].
An orthogonal approach to the selection of multiple template structures for sampling protein conformations is the generation of template protein ensembles by using MD simulations. For lactose permease (pdb code: 1pv6), a representative of the Major Facilitator Superfamily which to a high degree relates to uptake transporters by secondary structure, both all-atom [42,91] and coarse-grained [50] simulations have been performed for elucidating ligand binding and even transport events. However, direct usage of MD simulations to mimic the translocation process of template-based comparative models is only rarely used since the majority of available template structures cover only transmembrane regions. For example, available templates for OATP1B1 cover b400 amino acids out of 691 in total. Running MD simulations with an incomplete target structure could lead to artifacts. Furthermore, the limited time scale for all-atom simulations (hundreds of nanoseconds to a few microseconds), as well as the choice of including or omitting the lipid Fig. 3. Proposed computational workflow for studying ligand interactions with ADMET-relevant SLC transporters. Results from in vitro validation can be inputted to the stages of data curation and data analysis and subsequently be used for a new iteration of modeling. bilayer into the simulated system, can heavily affect the correctness and interpretability of the MD simulations. Contemporary simulation techniques, such as stochastic Monte Carlo, or MD simulations with enhanced sampling (e.g. by applying replica-exchange methods) can be used for structural refinement of extra-and intracellular domains.
To date, enhanced sampling techniques have been used only once for building a complete structural model for any of the discussed uptake transporters. As shown in a MD simulation study on OAT1 [112], modeling of the complete OAT1 transporter structure was divided into two stages: First, only the transmembrane region was modeled based on a template with high secondary structure similarity (glycerol-3-phosphate transporter; pdb code: 1pw4). At the second step, the extracellular domain was iteratively sampled by simulated annealing, while the transmembrane region was kept restrained. Results of this study indicate that the structural refinement by applying enhanced sampling methods could significantly improve existing structural models for uptake transporters which in turn would enhance the understanding of functional aspects of the transport mechanisms.
Only recently, reduced representation methods, such as normal mode analysis by applying elastic network models, can be used to overcome shortcomings arising from the high computational demands of conventional MD simulations ( [142,143]). For example, elastic network models have been used for the structurally-related fucose transporter (pdb code: 3o7q) to study its molecular basis for allosteric modulations [18]. Furthermore, by building comparative models in different functional states (inward-and outward-open conformation), elastic network models were capable to reproduce the whole translocation pathway of this transporter.
For multiple ensemble docking, normal mode simulations have been shown to be particularly useful to e.g. detect a biologically relevant conformation of dopamine D3 receptor, which has subsequently been prospectively validated by the existing dopamine D3 crystal structure [15]. One might argue that the generation of multiple conformations for a template structure might inherently bias the construction of structural models, since they can possess its internal dynamics which probably cannot be completely captured by the template. However, as shown and discussed in case of the human MATE1 homology model [130], the overall conformational stability between the template and the derived MATE1 homology models remained unchanged.
As a conclusion, the use of normal mode simulations in structurebased modeling studies can potentially improve conformational sampling when modeling uptake transporters. This in turn can lead to more accurate docking poses with the aim to better understand ligand-protein binding events and potentially selectivity switches.

Summary, Conclusions & Future Perspectives
ADMET-relatedSLC transporters are proteins of emerging interest in the framework of preclinical drug design. As demonstrated herein -by collecting available ligand and protein information from the open domaindata sparseness resulted in quite limited understanding of these transporter to date. Other factors complicating effective exploration of this class of proteins is their promiscuous nature, with potentially multiple binding sites, as well as overlapping substrate-and inhibitor profiles.
As demonstrated by discussed examples of molecular modeling and data analysis herein, new emerging technologies are on the rise also for these targets being particularly hard to unlock. Especially, data integration techniques and data analysis can lead to useful hypothesis about interesting SAR series at the beginning of an in silico study. Further, combining LB and SB methods seems to be an effective strategy, especially when it comes to selectivity profiling (like in the case of PCM modeling), or the exploration of knowledge about multiple binding sites (like in the case of combinatorial pharmacophores). In general, inclusion of in vitro experiments is a must especially for SB methods, e.g. to test the established binding mode hypotheses. In return, those in vitro measurements will lead to an increase in data points for a particular target, which can further be explored by statistical methods (such as machine learning approaches). For SB approaches, it would be interesting to include more systematically conformational sampling of protein conformations and multiple template structures into the comparative modeling step. It will be interesting to then compare results to those from docking into a single static template.
We are proposing a general workflow for in silico modeling of clinically relevant SLC transporters (see Fig. 3) which makes use of all available molecular modeling approaches and combines them with timely data science approaches (as far as the available data allows the different methods).
With the aim to develop safer medicine, it is of extraordinary importance to increase our understanding of the molecular basis of respective transporter-ligand interactions on an individual level (transporter by transporter), but also from a global point of view (how they act in concert). During evolution, transporters were optimized to help metabolizing chemical matter that we were exposed to. Obviously, they found a way to transport a wide variety of chemically distinct compounds. Understanding the interplay of clinically relevant transporters in transporting chemical matter, as well as the mechanisms underlying transporter selectivity, can help to unravel potential drug-drug or drug-food interactions which will finally lead to safer drugs in the future.