In silico optimization of pharmacokinetic properties and receptor binding affinity simultaneously: a ‘parallel progression approach to drug design’ applied to β-blockers

The present work exploits the potential of in silico approaches for minimizing attrition of leads in the later stages of drug development. We propose a theoretical approach, wherein ‘parallel’ information is generated to simultaneously optimize the pharmacokinetics (PK) and pharmacodynamics (PD) of lead candidates. β-blockers, though in use for many years, have suboptimal PKs; hence are an ideal test series for the ‘parallel progression approach’. This approach utilizes molecular modeling tools viz. hologram quantitative structure activity relationships, homology modeling, docking, predictive metabolism, and toxicity models. Validated models have been developed for PK parameters such as volume of distribution (log Vd) and clearance (log Cl), which together influence the half-life (t1/2) of a drug. Simultaneously, models for PD in terms of inhibition constant pKi have been developed. Thus, PK and PD properties of β-blockers were concurrently analyzed and after iterative cycling, modifications were proposed that lead to compounds with optimized PK and PD. We report some of the resultant re-engineered β-blockers with improved half-lives and pKi values comparable with marketed β-blockers. These were further analyzed by the docking studies to evaluate their binding poses. Finally, metabolic and toxicological assessment of these molecules was done through in silico methods. The strategy proposed herein has potential universal applicability, and can be used in any drug discovery scenario; provided that the data used is consistent in terms of experimental conditions, endpoints, and methods employed. Thus the ‘parallel progression approach’ helps to simultaneously fine-tune various properties of the drug and would be an invaluable tool during the drug development process.


Introduction
The pharmacological properties of a drug depend not only on its biological affinity at the receptor sites, but also on its ADMET profile. The key to develop a successful drug is to optimize its pharmaceutical attributes. A high potential exists for the use of in silico approaches to minimize the attrition rate during the drug development process. Hence, we have applied a fusion of methodologies to yield a 'parallel progression approach' to simultaneously optimize the pharmacokinetics (PK) and pharmacodynamics (PD) using structural information, which limits expending valuable resources and time while awaiting the outcome of iterative experimental evaluations. The advantage of this approach is that the modified molecules can be further evaluated for PK and PD, and during each progressive cycle further modifications can be made for lead optimization. Thus, the re-engineered structures can be taken further through the drug development process.
β-adrenergic blockers show a broad range of therapeutic utility, being used to treat hypertension and a variety of associated cardiovascular diseases (Frishman & Saunders, 2011). β-blockers were introduced at a time when the role of PK in the drug design process was poorly understood. As a consequence, most β-blockers introduced into therapy do not have desirable PK profiles (Obach, Lombardo, & Waters, 2008). Thus, this class of drugs forms an ideal test case for the 'parallel progression approach' for simultaneous optimization of their PK parameters and β 1 -blocking activity. Since there is a paucity of enantioselective PK data of β-blockers in humans, and also because most of the β-blockers used in therapeutics are racemates, the reported PK parameters for racemates were used in this study. It is pertinent to mention that although different β-receptor subtypes and associated β-blockers with receptor selectivities are known, we focused on drug design for molecules interacting with the β 1 -receptor for evaluation of the in silico strategy. The approach can indeed be applied to other classes of molecules and their associated receptors, and extension of the methodology with other β-receptor subtypes will be reported in future.
The chemometric tool used for this purpose is hologram quantitative structure activity relationships (HQSAR) (Heritage & Lowis, 1999;Lowis, 1997). This is an ideal tool for the execution of the parallel progression approach because it has the distinct advantage of requiring only 2D input structures for generation of a hologram (an extension of the 2D fingerprint), thus avoiding the need for molecular alignment and conformation specification inherent in the 3D methods like CoMFA and CoMSIA (Doddareddy et al., 2004). This is especially important because the experimental data used in this study is based on racemates. Thus, assigning correct 3D conformations to the molecules would not be possible for 3D-QSAR. Moreover, the PLS regression method in HQSAR is particularly suited when the matrix of predictors has more variables than observations and when there is multicolinearity among descriptor values which is the case here. The principle of parallel progression requires a chemometric tool that can model both PK, as well as, pK i parameters efficiently, so as to allow iterative cycles of design and prediction, quickly and accurately. HQSAR satisfies all these requirements since it is fast, reproducible, and does not involve complex descriptor selection process while still encoding all possible molecular information. Most importantly, HQSAR helps to identify fragments and maximal common structure (the template common to all training set molecules) that make important contributions to activity/kinetics through color-coded contribution maps which accompany all model outputs. This dual ability to suggest molecular fragments that are important contributors to the property, as well as, to predict the property for untested molecules makes HQSAR the perfect tool for this study.
The molecules with proven β-blocking activity were structurally modified using the parallel progression approach, resulting in re-engineered β-blockers that exhibit optimum PK profile in terms of half-life. Although, these new molecules were predicted to have binding affinity (pK i ) comparable with that of the marketed β-blockers, we further conducted docking studies in order to evaluate the binding poses and verify that the re-engineered molecules bind in a manner similar to the established β-blockers. However, the absence of a crystal structure for the human β 1 -adrenergic receptor necessitated the building of a homology model in order to carry out the docking studies.
As mentioned before, the objective of the study is to evaluate whether a late stage attrition of lead compounds during drug development can be minimized. Hence, two important aspects which have been prominently implicated in the late stage failures of drugs were also studied, i.e. metabolism and toxicity. In order to get an understanding of the metabolic profile, the re-engineered molecules were subjected to in silico metabolic assessment by performing induced fit docking (IFD) on cytochrome 2D6, which has been reported as the primary metabolizing isoform for the aryloxypropanolamine class of drugs (Coleman, 2010;Hlavica, 2006;Vaz, Nayeem, Santone, Chandrasena, & Gavai, 2005). Understanding of druginduced toxicity is invaluable in contributing to early compound prioritization and selection (Fielden & Kolaja, 2008). Several different kinds of in silico methods in toxicology have been developed and applied in the pharmaceutical industry (Raunio, 2011). The toxicity profiles of the re-engineered molecules were evaluated with the aid of in silico predictive toxicological models.
Thus, the current work endeavors to develop an in silico 'parallel progression approach' as a tool which will enable simultaneous estimation of all attributes of various molecules in the initial stages of drug discovery. This, in turn, would help in the selection of leads with optimal properties that would have much better prospects of evolving into successful drug candidates.

Background of in silico approaches
The wide therapeutic utility of β-blockers and the need for candidate molecules with better PK profiles has resulted in a number of approaches in the past aimed at developing predictive QSPkR (quantitative structure PK relationship) models for this class of drugs. The methods that were employed in these studies include univariate and multivariate correlations between important PK parameters such as V d and Cl, with simple physicochemical properties such as lipophilicity and ionization constant (Hinderling, Schmidlin, & Seydel, 1984). Another study involved construction of an artificial neural network (ANN) system based on a data-set of ten β-blockers, for predicting the PK parameters from the octanol/water partition coefficient, pK a or fraction bound to plasma proteins (Gobburu & Shelver, 1995). Although, neural networks are powerful tools and develop complex and robust relationships with good predictive ability, in the absence of a direct equation, these do not lend themselves easily for identification and fine tuning of important structural features. Some studies which report methods to predict PKs of β-blockers in the data-set are as mentioned (i) a study employed support vector regression method and a multiple linear regression method to develop models for predicting systemic clearance and steady state V d for data-sets including some β-blockers (Gombar & Hall, 2013); (ii) a study reported general regression neural network for modeling blood-brain barrier penetration, binding to human serum albumin, and milk-plasma distribution for data-sets including some of the β-blockers (Yap, Li, & Chen, 2006); (iii) ANNs have been applied in many QSPkR studies (Talevi et al., 2011;Turner, Maddalena, & Cutler, 2004;Wessel, Jurs, Tolan, & Muskal, 1998). Most of the QSPkR reports have shown that lipophilic/hydrophilic and electronic properties are more important than steric characteristics in determining the PK behavior of compounds (Kubinyi, 1997;Mager & Jusko, 2002;Testa, Crivori, Reist, & Carrupt, 2000;Mayer & van de Waterbeemd, 1985). Other studies have reported lipophilicity, solubility descriptors, polarity, and molecular size descriptors to be important parameters (del Amo et al., 2013;Karalis, Tsantili-Kakoulidou, & Macheras, 2002).
To the best of our knowledge, there have been no reports on in silico methodology for β-blockers that affords a comprehensive preview of an array of parameters vital for the development of novel drug molecules. Thus, the present work aims at fulfilling the need for an approach that encompasses simultaneous optimization of various aspects of drug design such as PK (V d , Cl), PD (pK i ), metabolism, and toxicity-related attributes in an in silico environment using β-blockers as a test case. The major benefits of this method is that being an in silico approach, it will help in minimizing the time and resource-intensive experimentation on a large number of molecules thus enabling synthesis of selected molecules. This method helps to understand the myriad factors affecting PK and PD, and also aids in quick identification of fragments in the drug structure responsible for specific effects. Most importantly it enables the application of the understanding gained, by allowing incorporation of changes in the structure and finally the modified molecule's PK parameters and β 1 -blocking activity can then be predicted easily by the same approach. Thus, the in silico parallel progression approach to drug design lends itself as a flexible tool to tune the various aspects of drug design simultaneously and ensure a higher probability of creating a successful drug candidate.

Data-sets and computational details
For the present work two data-sets, for modeling the PK and receptor binding affinity (pK i ) separately, were compiled from literature; the molecules therein and their corresponding experimental PK parameters and pharmacodynamic attributes are listed in the supporting information Tables S1-S3. When data was obtained from multiple reports, the average value was used. Data-set-1 (Table S1) encompasses the PK parameters for β-blockers. It consists of 32 diverse β-blockers with human experimental PK parameters i.e. volume of distribution (V d ) and clearance (Cl).
Data-set-2 (Table S2-S3) covers inhibition constants for β-adrenergic blockers expressed in terms of pK i values. Since many variables influence K i values such as species, cell lines, and functional or binding assay, the molecules were segregated based on the radiolabeled substrate i.e. [ 125 I]-iodocyanopindolol (ICYP) or (4-(3-tertiarybutylamino-2-hydroxypropoxy)benzimidazol-2-one) (CGP12177) used in the binding studies. Therefore, data-set-2 was further segregated into two classes, the first class comprising 43 compounds (PD data-set-I) is listed in Table S2 and the second class of 53 compounds (PD data-set-II) is listed in Table S3. Separate models were developed for the two classes.
All computational studies were carried out on a high performance computer cluster with 16 GB physical memory on each node running under the CentOS platform. The computations were executed with the molecular modeling software -Sybyl X Suite; HQSAR v1.1 (Tripos Inc., USA), Maestro, and ligprep modules (Schrödinger Inc., 2012, USA).
All molecules were prepared using the ligprep module and all possible states were generated at pH 7.0 ± 2.0. The β-blockers were found to be ionized at this pH. Of the various stereoisomers generated, the S-isomer was retained for aryloxypropanolamines and the R-isomer for the arylethanolamines.

Hologram quantitative structure activity relationships
HQSAR is the central chemometric tool used for modeling the PK and PD properties. It is essentially a QSAR technique that has the advantage of requiring only 2D structures as the input. It is based on the premise that the structure of a molecule is the determinant of all molecular properties. HQSAR generates an extended form of fingerprint, known as a molecular hologram that contains all possible molecular fragments within a molecule, including overlapping fragments, and maintains a count of the number of times each unique fragment occurs. Fragmentation of the molecules is done on the basis of various fragment distinction parameters selected. The 'atom distinction parameter' (A) distinguishes fragments on the basis of differences in their elemental types. The 'bond distinction parameter' (B) uses the differences in the bond types to differentiate fragments. Information about the hybridization state of the atoms in the fragments is incorporated in the 'atomic connection parameter' (C). The 'hydrogen inclusion/exclusion parameter' (H) differentiates fragments based on whether hydrogen atoms are included or excluded in the sybyl line notation strings generated. The 'donor and acceptor' (D) option initiates the search for donor and acceptor atoms in the fragments generated as defined. The 'chirality' option differentiates fragments on the basis of atom chirality and bond stereochemistry (Lowis, 1997). For each model, contribution maps can be generated. These are color-coded structural representations of the molecules. The color of each atom depicts its contribution to the overall molecular property. Atoms that negatively contribute to the property are represented by colors at the red end of the spectrum (red, red-orange, and orange) while those atoms that show favorable or positive contributions are represented by colors at the green end of the spectrum (yellow, green-blue, and green). Atoms with intermediate contributions are colored white. HQSAR models were generated for PK i.e. V d , Cl, and for PD properties i.e. pK i of the β-blockers using various combinations of the fragment distinction parameters and hologram lengths.

Modeling PKs and receptor binding affinity
During drug development process, the PK and PD aspects of a molecule are the critical determinants of its potential to develop into an effective drug candidate. The ability to identify molecular determinants which affect the PK and pK i profiles in the initial stages of drug design would be highly desirable. Towards this end, in the present work, fragment based in silico QSPR models have been developed to predict PKs and receptor binding attributes of β-blockers.

PK model for volume of distribution (V d )
A combination of different hologram lengths and various fragment distinction parameters was used to generate the HQSAR models. Nineteen β-blockers comprised the training set and four were in the validation set. The model with the best statistical parameters is reported here, where the fragment distinction was done on the basis of atoms/bonds/connections/donor or acceptor (A/B/C/D) parameters with fragment length 5-9 (F size ), hologram length 53 (HL), and 5 components.

PK model for clearance (Cl)
The training set comprised of sixteen molecules, and four molecules were in the validation set. The model with the best statistical parameters was obtained using atoms/bonds/donor or acceptor (A/B/D) as the fragment distinction parameters, F size 4-6, HL 83, and 4 components.

PD model for inhibition constant (pK i )
The fragment distinction parameters that gave the best result for PK data-sets-I and II were atoms/bonds/connections (A/B/C). The training set comprised of 27 molecules and eight molecules were in the validation set for PK data-set-I, wherein the best model was developed with F size 4-7, HL 401, and 6 components. Thirty molecules were included in the training set and eight molecules in the validation set for PD data-set-II, and the best model developed with F size 7-10, HL 61, and 6 components.

Receptor binding studies
The HQSAR methodology yielded predictive models for the various PK and receptor binding affinity (PD) parameters for β-blockers, and these were used to re-engineer the structure of the molecules. In order to gain a detailed understanding of the binding mode and interactions of the established and modified (re-engineered) β-blockers to the β 1 -adrenergic receptor, docking studies were carried out.
The crystal structure of the human β 1 -adrenergic receptor has not yet been reported, though the crystal structure of the turkey β 1 -adrenergic receptor is available (2VT4) (Warne et al., 2008). Hence, homology modeling was done to generate an in silico model of the human β 1 -adrenergic receptor, which in turn was used to study the binding poses of known β-blockers as well as the re-engineered β-blockers.

Homology modeling
A number of homology models were generated by submitting the primary sequence of the human β 1 -adrenergic receptor with various web-based servers such as ITASSER (Roy, Kucukural, & Zhang, 2010;Zhang, 2008), GPCR-SSFE (Worth, Kreuchwig, Kleinau, & Krause, 2011), GPCRDB, M4T-cluster (Fernandez-Fuentes, Madrid-Aliste, Rai, Fajardo, & Fiser, 2007), and GPCR ModSim (Rodríguez, Bello, & Gutiérrez-de-Terán, 2012). The models were built on templates with high sequence identity and similarity with the human β 1 -adrenergic receptor. The best model was selected on the basis of various parameters such as the Ramachandran plot, Maestro align score, verify 3D, and DOPE score. Since the β-blocker binding pocket is located in the transmembrane region of the receptor, the selected model was truncated at the N-and C-terminus to reduce computational time. The structure was then refined by molecular dynamics (MDs) simulation.

MDs simulation
Since the β 1 -adrenergic receptor is a G-Protein coupled receptor, it was simulated in a DPPC (dipalmitoylphosphatidylcholine) bilayer membrane. The coordinates for the membrane were obtained using the PPM server (Lomize, Pogozheva, Joo, Mosberg, & Lomize, 2012). The protein structure was minimized using the OPLS 2005 force field till a gradient of 0.1 kcal/mole/Å. MD simulations were carried out using Desmond v3.1 (Shaw, 2005) in the Schrödinger suite 2012. The OPLS-2005 force field was used and the system was solvated using the TIP3P water model (Mark & Nilsson, 2001). Appropriate amounts of Na + and Cl − ions were added to maintain the salt concentration of the system at 0.15 M, to maintain physiologically relevant concentration (2012). This solvated complex was relaxed with restraints on the solute and subsequently equilibrated at 325 K and 1 atm pressure. Later, the equilibrated system was simulated for 50 ns under NPT ensemble to study the time dependent fluctuations and conformational changes of the receptor model. The structural states and the corresponding energies were sampled from the trajectories every 50 ps.

Docking
Molecular docking studies were carried out to predict and compare the 'most probable' bound association of the existing and re-engineered β-blockers with the homology model of the human β 1 -adrenergic receptor. These studies were performed using GLIDE 5.8 in Schrödinger Suite 2012 running on an Intel Xeon processor-based workstation with the CentOS enterprise Linux 5.5 (ROCS cluster 6.1) OS. Initially, the homology model of the β 1 -adrenergic receptor was prepared using the protein preparation wizard, wherein charges, atom types, bond orders, and formal charges were correctly assigned, and the protein terminus was capped with ACE and NMA residues. The protein's hydrogen bond network was optimized using a systematic, cluster-based approach and restrained minimization was carried out using the OPLS 2005 force field that allowed hydrogen atoms to be freely minimized, while allowing sufficient heavy atom movement to relax strained bonds, angles, and clashes such that the heavy atom positions converged to an RMSD of 0.1 Å.
In order to setup a validated docking protocol, native ligands of the turkey β 1 -adrenergic receptor such as cyanopindolol, carazolol, and carvedilol were docked to the homology model. The receptor grid of size 20 Å was generated by taking the centroid of the native ligand cyanopindolol as the grid center. The conformations were selected from an exhaustive enumeration of the minima in the ligand torsion-angle space and prescreened over the entire phase space available to the ligand to locate promising ligand poses, and minimized in the field of the receptor using the OPLS 2005 force field in conjunction with a distance-dependent dielectric model. The lowest energy poses obtained in this fashion were subjected to a Monte Carlo procedure that examines nearby torsional minima (Friesner et al., 2004). The docking protocol was run and the poses were ranked on the basis of the Glide 'standard precision' docking score. The ligand binding interactions were analyzed and compared with those seen in the crystal structure. The re-engineered β-blockers with optimal PK and receptor binding affinity (PD) properties were docked using the validated protocol.

In silico metabolism studies of re-engineered β-blockers
Cytochrome P450 enzymes play an integral role in the oxidative metabolism of drugs including β-blockers. The ability to predict the reactivity and accessibility to CYP isoform and the possible sites of metabolism is necessary to understand the metabolic stability of re-engineered molecules. This was studied for various marketed β-blockers and the re-engineered ones by 'induced fit docking' on CYP 2D6 isoform using the 'P450 Site of Metabolism' module of Schrödinger Suite v9.3. For a given atom in a molecule to be a significant site of metabolism by a CYP450 isoform, it must have some degree of 'reactivity' in the absence of the enzyme and also be accessible to the reactive heme iron center. To address both these requirements, the 'P450 Site of Metabolism' workflow, which combines IFD for the determination of accessibility to the reactive center, with a rule-based approach to intrinsic reactivity, was used.

In silico toxicity studies
For this study, ADMET predictor v6.0 implemented in the Simulation Plus, Inc., software was used to predict the toxicity profiles of the re-engineered molecules and these were compared with existing β-blockers. ADMET predictor utilizes QSAR-based expert systems that are mainly used in early drug discovery to predict toxicological endpoints such as carcinogenicity, teratogenicity, mutagenicity, immunotoxicity, neurotoxicity, developmental toxicity, respiratory sensitization, and skin irritation etc.

Modeling PKs and receptor binding affinity
The statistics for the best models are summarized in Table 1. Beside the standard parameters q 2 , r 2 pred , two additional parameters r 2 m and R 2 p were also calculated (Pratim Roy, Paul, Mitra, & Roy, 2009).

PK Models for V d and Cl
The contribution of atoms toward defining the V d is shown as contribution maps in Figure 1. Analysis reveals that the presence of a sulfonamido substituent negatively contributes to the V d , whereas presence of a cyclopropylmethoxy ethyl side chain does not contribute significantly towards V d .
The atomic contributions toward defining clearance can be visualized in Figure 2. In contrast to the V d model; the atomic contributions defining Cl illustrate that in most cases the central aromatic moiety makes a positive contribution, as also the presence of methyl substituents on the aromatic ring.
All models have been validated using molecules with known experimental values and which have not been a part of the training set for the construction of the models. The results are depicted in supplementary Table S4.

Modifications on existing β-blockers
Since the existing β-blockers have ample scope for improvement in their PK profile, structural modifications have been made based on the contribution maps obtained for the V d and Cl models. Careful analysis of the contribution maps revealed certain structural handles that could be used to increase/decrease V d and/or Cl. There is no optimal value for either V d or Cl; however, these together influence the half-life (t 1/2 ) of a molecule. Formulation and patient compliance aspects, suggest that a half-life between 6 and 10 h is considered to be ideal. Hence, the strategy was to incorporate structural modifications to alter V d and Cl and thereby influence the elimination rate constant (k), where k = Cl/V d (the halflife is related to k as t 1/2 = 0.693/k), so as to fine-tune the half-life of the re-engineered molecules between 6 and 10 h.
Cetamolol was selected as the candidate for structural modification as it has median values for both V d and Cl and lends itself easily to modification. The contribution maps for V d for most molecules showcase fragments that make a negative contribution hence the strategy for increasing V d for these molecules was to eliminate such substituents. V d was increased by incorporating structural handles such as a cyclopropylmethoxy ethyl side chain on the aromatic ring or N-methylaminocarbonyl side chain or changing the position of the side chains such that they are in a 1-, 4-disposition to each other. Few structural handles viz. addition of methyl groups on the ring, replacing the N-methylaminocarbonylmethoxy side chain by acetoxy or 2-(morpholinocarbonylamino) ethoxy, introducing an additional butanamide side chain and replacement of the benzene ring with an indole nucleus lead to an increase in clearance. The structural modification suggested by the V d contribution maps may alter the Cl profiles and it is necessary that 'structural modification,' as well as, 'V d and Cl predictions' are done simultaneously and analyzed after each cycle of modification. Only those structural changes were retained, that leads to re-engineered molecules with halflives between 6 and 10 h, as depicted in Table 2. Though numerous structural modifications are theoretically possible, the modifications proposed here result from fragments which have been extracted from the HQSAR contribution maps. The metabolic susceptibility, toxicity,    and synthetic feasibility were important determinants in their selection.

Modeling pK i
HQSAR models for pK i of molecules in the two datasets (PD data-set-I and II) were developed. The statistics for the best pK i models are also showcased in Table 1.
The atom-wise contributions toward the inhibition constant (pK i ) for the various molecules are depicted in Figures 3 and 4 for molecules in PD data-sets I and II, respectively. The analysis of models thus developed reveals that the 'aryloxypropanololamine' fragment reported to be vital for the β-blocking activity, was identified as the maximal common substructure (colored cyan). Modifications in the   aryl portion such as changing the phenyl ring into a naphthyl, addition of methyl groups to the aryl system showed a positive contribution to pK i .
3.1.4. Parallel progression approach for PK and receptor binding affinity (PD) modeling The next step was to simultaneously optimize the PK and receptor binding affinity (PD) profiles. A concurrent analysis of the HQSAR contribution maps for the PK and pK i (PD) models of the existing β-blockers was carried out to identify fragments that could improve β 1 -adrenergic blocking activity while simultaneously optimizing the PK profile of the molecules (t 1/2 ). These features were incorporated to yield re-engineered molecules with the desired product profile. The pK i , V d , and Cl parameters of the re-engineered molecules were predicted using the previously built and validated HQSAR models. Progressive modifications based on iterative cycles of PK and pK i (PD) models led to modified molecules with the optimal PK and pK i values. A few of these re-engineered molecules are depicted in Table 3. The PK and receptor binding affinity studies were further fortified by the following studies.

Homology modeling and MDs simulation
Multiple sequence alignments of the human β 1 -adrenergic receptor using Clustal Omega (Sievers et al., 2011) revealed that the turkey β 1 -adrenergic receptor (2VT4) with 75% identity could be used as a suitable template for homology modeling. Various web-based servers were used to develop homology models for the human β 1 -adrenergic receptor. After truncation at the N-and C-terminals, the models were aligned with the crystal structure of 2VT4. The structural stability of the models was assessed by comparison of the predicted secondary structure with the crystal structure, and on the basis of the RMSD with reference to the template. The homology model generated using the GPCR ModSim server is based on 2VT4 as the template. This was selected for further studies as it showed the least RMSD with a maestro align score (all atoms) of 0.153, align score (helix only) of 0.004 on alignment with the crystal structure (2VT4); DOPE (discrete optimized protein energy) score (Shen & Sali, 2006) of -40309.24, profile 3D verify score (Velikanov, Yan, Badretdinov, & Szalma, 2001) as 75.54 and only 2% of the residues in the disallowed region in the Ramachandran plot (Hooft, Sander, & Vriend, 1997). It was observed that the residues involved in ligand binding are conserved. A detailed analysis of the secondary structure composition was performed using the PDBSum tool (Laskowski, Chistyakov, & Thornton, 2005) and it was verified that the homology model has similar structural motifs and overlapping transmembrane regions as the crystal structure (2VT4). MD simulation for 50 ns was performed on this homology model embedded in DPPC, a membrane mimetic, to allow the protein to relax from its initial conformation. The Simulation Event Analysis of the trajectory using visual molecular dynamics (Humphrey, Dalke, & Schulten, 1996) was used to identify the structure to be used for the docking studies. The last frame in the MD simulation was selected because the RMSD plot showed minimal atomistic fluctuations of the backbone atoms. Scrutiny of the Ramachandran plot indicated that no residues were in the disallowed region as opposed to 2% prior to MD simulation. Further validation of the homology model was done by analyzing various parameters like G-factor (for stereochemical properties such as torsion angles and covalent geometry), bond lengths, and bond angles (analysis of the main-chain bond lengths and bond angles of the protein structure), planar groups, residue properties (bad contacts) which were calculated using PROCHECK suite of programs (Laskowski, MacArthur, Moss, & Thornton, 1993). Furthermore, ANOLEA (Atomic Non-Local Environment Assessment) (Melo & Feytmans, 1998) sever was used for calculating and evaluating the 'Non-Local Environment' (NLE) of each heavy atom in the molecule. The results are summarized in Table 4. Alignment of the binding site of the turkey β 1 -adrenergic receptor with the human β 1 -receptor homology model revealed that the orientation of Asn363 was different in the two cases, hence side chain refinement was done using Discovery studio (Studio, 2010). A rotamer search showed that the Ponder and Richards rotamer 2 for Asn was important to mimic the ligand receptor interactions. The homology model is depicted in Figure 5.

Docking
The important binding interactions observed in the crystal structures of ligands such as cyanopindolol, carazolol, and carvedilol bound to the β 1 -receptor include hydrogen bonds between the -OH group of the ligand and the carbonyl oxygen of Asp138 and between the oxygen atom of ligand's -OH group and the -CONH 2 of Asn363. The terminal amino group of the ligand is hydrogen bonded to the carbonyl oxygen of Asp138 and/or Asn363. Also, a π-stack is observed between the central aromatic moiety of the ligand and Phe341. These interactions were reproduced during the validation studies of the docking protocol. The validated protocol was then used to study the binding of the re-engineered β-blockers with the homology model of the β 1 -receptor. The results indicate that the binding modes of the re-engineered molecules are similar to that observed in the crystal structures of ligands bound to the β 1 -receptor, as seen in Figure 6.

In silico metabolic studies
A metabolic reaction can proceed at a highly reactive site of a molecule, only if it is properly oriented within the CYP enzyme to interact with the reaction center. Thus, intrinsic reactivity and accessibility are the decisive aspects for a metabolic reaction to take place. Hence, metabolism of the re-engineered β-blockers was studied by IFD with CYP2D6 and the results were analyzed on the basis of 'Fe accessibility,' 'intrinsic reactivity,' and 'overall SoM (Site of Metabolism)'. Details of these parameters, the metabolic maps of an established β-blocker as a reference and the re-engineered β-blockers are illustrated in Figure 7. The metabolic map for cetamolol shows high probability for N-demethylation as the N-methyl group shows both high Fe accessibility, as well as, good intrinsic reactivity. Also the maps indicate probability of aromatic hydroxylation due to their proximity to Fe however, their intrinsic reactivity is low. Similarly, though O-dealkylation and oxidative deamination to aryloxycarboxylic acid is possible as intrinsic reactivity is high, however Fe accessibility is quite low indicating that these metabolic pathways will yield minor metabolites. These predictions are in line with the experimental data available for cetamolol metabolism. A perusal of the predicted metabolic maps for EPB-4 and EPB-56 shows metabolic profiles similar to Cetamolol. One of the notable differences is that the N-isopropyl fragment of EPB-56 lacks Fe accessibility Figure 6. Comparison of the key interactions of established β-blocker. Carvedilol (A) and re-engineered molecules EPB 56(B), EPB 4(C) [depicted by the tube representation], with the homology model of the β 1 adrenergic receptor. Residues involved in H-bonding are shown in green and those involved in π-π stacking in purple. The hydrogen bonds are depicted in red. and hence, the N-dealkylation pathway will not be favored. However, the methyl cyclopropyl side chain shows high Fe accessibility and suggests probability of O-dealkylation.
Overall, it can be concluded that in silico the re-engineered molecules show intrinsic reactivity and accessibility for the CYP 450 system comparable to the known drugs and thus we infer that they would not be more susceptible to metabolism than their marketed counterparts.

In silico toxicity studies
The increased understanding of biological systems and the tremendous advances in computational modeling techniques has made it possible to develop in silico approaches that aid in early identification of the potential hazards of NCEs (new chemical entities), thus leading to the reduction of development risk and cost in the early phase of drug discovery. This study reveals that the re-engineered β-blockers are predicted to be nontoxic for the endocrine receptor and the hERG K + channel (which indicates the potential of a compound for cardiac toxicity). Although these molecules are predicted to be toxic for the androgen receptors, their affinity level for these receptors (e.g. molecule EPB_4 = 0.0058) listed in Table S5 is in the same range as of the established β-blockers (e.g. cetamolol = 0.0064). The re-engineered β-blockers were also predicted to be non-sensitizers in a qualitative assessment of allergenic respiratory sensitization in rat and turn out to be nontoxic in a qualitative estimation of the ability to trigger the mutagenic chromosomal aberrations. They are also predicted to be nontoxic in the qualitative estimation of reproductive/developmental toxicity.
Thus, it is observed that the re-engineered β-blockers show better or comparable toxicity profile as that of the existing β-blockers. The results are given in the supplementary Table S5.

Conclusions
The current work is an attempt to develop an innovative theoretical approach that can simultaneously optimize both the PK and receptor binding (pharmacodynamic) attributes of a drug. The HQSAR tool was selected for this task. The β-blockers were selected as a test case. HQSAR models were developed for the important PK parameters i.e. volume of distribution (V d ) and clearance (Cl). In parallel, models were also developed for the inhibition constant (pK i ) that reflect the pharmacodynamic attributes. It should be emphasized here that though this study has been based on human PK and PD data; the approach itself can be used for any set of data (human or animal) that is consistent in terms of experimental conditions, endpoints, and methods employed. The PK models and in conjunction the pK i (PD) models were perused to modify the existing β-blockers to produce re-engineered molecules with optimal half-life and binding affinity to the receptor. Further it was established that the re-engineered molecules possess favorable receptor binding, metabolic, and toxicological profiles. Of course experimental data is needed to corroborate the results obtained through in silico methodology. It must be stated here that models were built from PK data that was obtained from human clinical trials. Any experimental validation of the designed molecules would involve human subjects necessitating approvals from government agencies, as this involves legal and ethical issues associated with human subjects.
Although, there are instances of multiparameter modeling in the literature, the novelty and the genuine intent of this work lies in the fact that the approach has not terminated by simply suggesting validated models, instead the research has been carried forward by analyzing the model predictions and utilizing these pointers to design re-engineered molecules. Further, a holistic measure of the metabolism and toxicity of these molecules has been made through a variety of well-established prediction tools. Thus, the current work tries to address most of the variables in the dynamic drug design process right at the outset with the convenience of in silico tools in order to boost the potential of the selected molecules to serve as good leads in terms of optimum PKs, pharmacodynamic, and toxicological attributes.
In conclusion, such a comprehensive in silico 'parallel progression approach' encompassing all aspects of drug design would help in making critical 'gono go' decisions in the early phases of drug development and thus improve the prospects for compounds in trials.

Supporting information
The supporting information includes Table S1: PK data for β-blockers, Table S2: pK i values of β-blockers determined using ICYP as the radiolabeled substrate, Table S3: pK i values of β-blockers determined with CGP12177 as the radiolabeled substrate, Table S4: Validation sets for pharmacokinetic and pharmacodynamic models, and Table S5: Results of toxicity studies performed with ADMET Predictor and references.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.org/10.1080/07391102.2015.1033646.

Disclosure statement
In accordance with Taylor & Francis policy, the authors declare that there is no conflict of interest regarding the research presented.