Deciphering ophthalmic adaptive inhibitors targeting RON4 of Toxoplasma gondii: An integrative in silico approach

Accepted Manuscript Deciphering ophthalmic adaptive inhibitors targeting RON4 of Toxoplasma gondii: An integrative in silico approach

Umashankar Vetrivel, Hemavathy Nagarajan PII: DOI: Reference:

S0024-3205(18)30639-8 doi:10.1016/j.lfs.2018.10.022 LFS 16002

To appear in:

Life Sciences

Received date: Revised date: Accepted date:

21 August 2018 9 October 2018 11 October 2018

Please cite this article as: Umashankar Vetrivel, Hemavathy Nagarajan , Deciphering ophthalmic adaptive inhibitors targeting RON4 of Toxoplasma gondii: An integrative in silico approach. Lfs (2018), doi:10.1016/j.lfs.2018.10.022

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Title Page Title: Deciphering ophthalmic adaptive inhibitors targeting RON4 of Toxoplasma gondii: an

Centre for Bioinformatics, Kamalnayan Bajaj Institute for Research in Vision and

CR

1

IP

Authors: Umashankar Vetrivel1*+ & Hemavathy Nagarajan1+

T

integrative in silico approach

Tamil Nadu, India. Telephone: +91-44-28271616

AN

E-mail: [email protected], *Corresponding Author

CE

PT

ED

M

Contributed equally

AC

+

US

Ophthalmology, Vision Research Foundation, Sankara Nethralaya, Chennai - 600 006,

ACCEPTED MANUSCRIPT Abstract Aims: Toxoplasma gondii (T. gondii) is an intracellular Apicomplexan parasite and a causative agent of toxoplasmosis in human. Its parasitic invasion leads to encephalitis, uveitis, chorioretinitis and congenital infection. Host invasion by T. gondii is mediated by Moving

T

junction (MJ) complex formed by secretory proteins such as micronemes and Rhoptry Neck

IP

proteins (RONs). Among these secretory proteins, RON4 shows major interactions with RON2

CR

and RON5 of MJ complex and is also found to facilitate invasion by interacting with the

US

tubulin of the host. Therefore, drug targeting of RON4 is considered to be an effective modality in inhibition of host invasion by T. gondii.

AN

Main methods: Hence, in this study, we have implemented High throughput virtual screening (HTVS) against predicted novel druggable cavity on RON4, in order to prioritize potential

M

inhibitors. This study also specifically focuses on identifying potential drugs that are ocular

ED

accommodative towards targeting ocular toxoplasmosis. Thus, the small molecules resulting

PT

from HTVS were subjected to stringent multi-level precision screening which involves PASS prediction, Density functional theory analysis, Binding Free Energy Calculations and stability of

CE

complex formation during molecular dynamics simulation to prioritize the potential hits. Key findings: Among the compounds from NCI chemical library, five (338683, 333739, 686585,

AC

159706, 405935) were found to surpass all the stringent filtration criteria. Significance: Based on comparatively analysis, it could be inferred that 333739 (benzyl 2-((2(((benzyloxy) carbonyl) amino)-3-hydroxybutanoyl) amino) propanoate) to be highly potential in comparison to other compounds and shall prove to be an efficient inhibitor of RON4 in ocular toxoplasmosis. Keywords: Moving Junction Complex, Toxoplasma gondii, RON4, virtual screening

ACCEPTED MANUSCRIPT 1. Introduction Toxoplasmosis is a chronic infective disease condition among humans. It majorly occurs as congenital infection and is also more prevalent in immuno-compromised individuals. Toxoplasma gondii is an obligate Apicomplexan intracellular parasite which causes

T

toxoplasmosis [1]. In case of congenital infection, it leads to acute ocular toxoplasmosis resulting

IP

in focal Chorioretinitis and other severities [2]. Many of the currently practiced are limited in

CR

efficacy and also confer toxic side effects [3–5]. Moreover, treatment of ocular toxoplasmosis demands more specific and sensitive drugs to combat the raising antimicrobial resistance,

US

thereby, posing a need for new drugs effectively modulating the pathogenesis. In T.gondii, host

AN

invasion is facilitated by Rhoptry neck (RON) and microneme proteins which anchors the parasite to the host membrane through the formation of moving junction (MJ) complex [6,7].

M

Targeting the proteins that form MJ complex has been proven to be effective modality for

ED

targeting the host invasion by T.gondii. The major MJ Complex proteins involved in host invasion are RON2, RON4, RON5 and RON8 [8,9].

PT

Among the RON proteins, only RON4 is found to commonly interact with RON2 and

CE

RON5 during the formation of MJ complex. Moreover, RON4 is also reported to interact with the host cellular β-tubulin [10,11], thereby, acting as a medium for localizing the parasite to the

AC

host. In addition, previous studies also highly signify that RON4 might function independent of apical membrane antigen1 (AMA1) [12] and other RON/ROP proteins to invade host [13]. RON4 being a cytosolic protein, is secreted by T.gondii on host cell surface during invasion. The secreted RON4 forms complexes with other RON proteins and aids in reorientation of the entire parasite. This leads to the close proximal assembly of these secretory proteins to get accumulated at apical region, thereby, establishing contact with the host [13].

ACCEPTED MANUSCRIPT Moreover, recent studies also infer that RON4 based drug targeting could be more potential in terms of specificity, as RON4 is non-homologous to human proteins [14]. Therefore, targeting RON4 on its druggable region would be an effective approach in order to disassemble the MJ formation, thereby, inhibiting host invasion. As RON4 lacks experimentally determined crystal

T

structure, a plausible 3D structure of RON4 was predicted in our previous study by

IP

implementing fold recognition and molecular dynamics simulation approaches. The predicted

CR

structure of RON4 (948 aa in length) comprised of 60 alpha helices, 58 beta turns, 15 gamma turns and 119 helix-helix interactions [14]. In our previous study, a potential inhibitor peptide

US

was proposed based on intermolecular interactions of host tubulin with RON4. However, there

AN

is a paucity of studies on identifying novel druggable regions in RON4 and targeting the same with small molecules. Molecular modelling tools and Computational aided drug design have

M

accelerated the drug discovery process. Effective and integrative usage of these tools shall lead

ED

to potential molecules for therapeutic applications [15,16]. Hence, this study was framed to

PT

predict potential RON4 inhibitors that are ocular accommodative by integrating High Throughput Virtual screening, Molecular dynamics simulation, quantum mechanics calculations,

AC

CE

Binding free energy calculation and machine learning approaches.

ACCEPTED MANUSCRIPT 2. Materials and methods 2.1 Binding site prediction The plausible 3D structure of RON4 predicted in our previous study [14] was used for mapping the novel druggable regions, wherein, the regions on RON4 interacting with AMA1 (C-

T

terminal) and -tubulin were excluded (727-860). Therefore, the predicted structure of RON4

IP

was subjected to binding site prediction using SiteMap module of Schrödinger suite [17,18]. The

CR

SiteMap module of Maestro operates similar to GRID algorithm of Good ford, wherein, the

US

calculations are performed at three stages. As a first step, a grid is set over the protein surface and the grid points are grouped into sets based on Van der Waals and distance-dependent

AN

electrostatic-interactions. Further, the identified sets of points were used to generate surface maps. Finally, the properties of the predicted sites based on surface maps were evaluated using

M

SiteScore and Druggability score. SiteScore is calculated as,

ED

SiteScore = 0.0733 sqrt (n) + 0.6688 e - 0.20 p

PT

Druggability Score is calculated as follows, Dscore = 0.094 sqrt (n) + 0.60 e - 0.324 p

CE

Where n denotes the number of site points, e is the enclosure score, and p is the hydrophilic score that has been set to a maximum of 1.0 based on the impact of hydrophilicity at the charged and

AC

highly polar sites. Apart from these scores, SiteMap also provides other properties like hydrophilic and hydrophobic nature, balance, donor/acceptor character, site volume etc in order to classify the better binding site. For an ideal drug binding site, Site Score greater than 1 is preferred. In case of druggablility score, the optimal druggable region will have a score nearer to 1. Thus, in this study, SiteMap analysis was performed to identity the top ranking druggable binding site and was further used for virtual screening studies.

ACCEPTED MANUSCRIPT 2.2 Virtual screening-RON4 2.2.1 Protein preparation and Receptor grid generation The minimum free energy conformation of RON4 [14] was pre-processed and optimized using Protein Preparation Wizard of Schrödinger (New York, NY, USA) [19]. The stereochemistry of

T

RON4 structure was optimised by assigning proper bond order, removing steric clashes. In

IP

addition, the optimal protonation states for histidine residues were also assigned and the chi

CR

rotation of aspargine, glutamine and histidine residues were also adjusted accordingly. Finally, energy minimization was performed with OPLS_2005 force field to obtain the optimal geometry.

US

Further, the top ranking druggable binding site on RON4 (receptor) predicted by SiteMap

AN

module was used for docking grid generation. Following which, the receptor grid was set with vanderwaals radii scaling factor set to 1.0 Ǻ and a partial cut off of 0.25 Ǻ and proceeded for the

M

virtual screening process.

ED

2.2.2 Ligand preparation:

To facilitate virtual screening of ligands from National Cancer Institute (NCI) repository

PT

(2,51,616 compounds) [20,21], the ligands were screened for the drug-likeness by implementing

CE

the Lipinski rule. The filtered drugs like molecules were further optimized by Ligprep module of Schrodinger suite [22] in order to fix the ring isomers, stereo isomers and tautomeric forms.

AC

Further, these structures were energy minimized using the OPLS 2005 force field. Subsequently, the ligands were also filtered for desirable physiological properties by applying Qikprop module of Schrodinger suite. In addition to this, ligands with reactive functional groups were also removed. Further, the pruned dataset was scaled with van der Waals radii scaling factor of 0.80 Ǻ with partial charge cut-off set to 0.15 Ǻ. Finally, these optimized ligand datasets were proceeded for virtual screening against the druggable site of RON4.

ACCEPTED MANUSCRIPT 2.2.3 Virtual screening workflow High throughput virtual screening (HTVS) and docking of optimized NCI compounds (2,51,616 compounds) against RON4 were performed using the Glide HTVS of the Schrodinger suite 2017 (Schrödinger, LLC) [23–25]. As a first step, High-throughput virtual screening (HTVS) was set

T

to generate 100 poses and top 10% of the compounds were prioritized based on the binding

IP

affinity. This was followed by standard precision (SP) docking, wherein, 100 poses per ligand

CR

(for the top 10% for HTVS) were generated and 10% of the compounds were subjected to extra precision docking (XP). Finally, 10% of the compounds obtained from extra precision docking

US

(XP) were analysed for its docking energy, docking score and binding free energy calculation.

M

were chosen for further in silico validation.

AN

Among the XP resultant compounds, the top 10 compounds with significant binding affinity

ED

2.3 Prediction of Activity Spectra of Substances (PASS) prediction analysis of hit compounds

PT

The top hit compounds resulting from XP mode were subjected to various filtration steps.

CE

Initially, these compounds were exposed to Prediction of Activity Spectra of Substances (PASS) analysis [26,27]. As PASS describes the intrinsic properties of compounds based on the

AC

biological activity spectrum (with 95% accuracy). The compounds with activity of treating ophthalmic ailments were further shortlisted and were subjected to molecular mechanics generalised born and surface area solvation (MMGBSA) calculations [28] to infer the binding free energy.

ACCEPTED MANUSCRIPT 2.4 Prime based Molecular Mechanics Genralised Born Surface Area (MMGBSA) scoring: MMGBSA scoring was performed for PASS qualified complexes, in order to calculate the binding free energies of the complex (with implicit solvent) by implementing Prime module of Schrodinger suite [29]. MMGBSA score is calculated in three dimensional space based on the

IP

T

following equation:

CR

△GGB = △EMM + △Gsolv + △GSA

Where △GGB represents the electrostatic solvation energy, △EMM denotes the difference

US

between the minimized energies of ligand-protein complex and the total energies of protein and ligand in free form. △GSA is the difference in the surface area energies for the protein and the

AN

ligand, whereas △Gsolv is the difference in the GSA of the ligand-receptor complex and the sum of

M

their solvation energies in the unbound state. The higher negative value of MMBGSA score

ED

indicates stronger binding conformation; therefore based on this, top ranking compounds were identified and subjected for MD simulations in order to further determine the stability of

PT

complex formation [30]. The intermolecular interactions of all the potential compounds with

CE

RON4 protein were validated by Ligand interaction module of Schrodinger suite (2017) and also

AC

using protein-ligand interaction profiler (PLIP) [31].

2.5 Density functional theory (DFT) analysis: DFT is a Quantum mechanics based approach which aids in understanding the molecular orbitals of drug like molecules, thereby inferring chemical reactivity and stability of the compounds. Thus, all the identified ophthalmic adaptive drugs in complex with RON4 were subjected to DFT analysis to determine molecular orbital energies. Among the molecular orbitals, HOMO refers to highest occupied molecular orbital and LUMO being the lowest

ACCEPTED MANUSCRIPT unoccupied molecular orbital (LUMO). HOMO energy defines the electron donor region of small molecules during complex formation, while LUMO energy denotes electron acceptor region. The difference in energy level between HOMO and LUMO is known as HOMO-LUMO energy gap (HLG), which signifies the electronic excitation energy that highly contributes for

T

computing the molecular reactivity and stability of the compounds. In this study, Jaguar module

IP

of Schrodinger suite [32] was used to optimize the geometry of the top 5 docked ligands based

CR

on Beckes three-parameter exchange potential and Lee-Yang-Parr correlation functional (B3LYP) theory with 6-31G** basis set. The basis set generally refers to the set of one-particle

US

(non-orthogonal) functions that acts as a base to build molecular orbitals. In the usage of 6-

AN

31G** basis set, a split-valence double-zeta plus polarization along with polarization for hydrogen atoms was utilized. Furthermore, HOMO- LUMO energy levels were calculated based

M

on surfaces such as molecular orbital, density, potential and atomic electrostatics potential

ED

charges (EPS). Finally, the reactivity and stability of the small molecules were inferred based on

PT

the minimal difference between HLG levels[33].

CE

2.6 Refinement and validation of the docked complexes using Molecular Dynamics (MD) Simulation

AC

The stability RON4 structure (Apo form) under MD simulation was confirmed in our previous study [14], hence it was not repeated here. Therefore, only RON4 in complex (Holo forms) with top five compounds were subjected to MD Simulation to sample the near native conformations, in order to validate the stability of complex formation. The simulations were performed using Gromacs 5.0.4 package [34,35] with OPLSAA (All-atom optimized potentials for liquid simulations) as force field. To start with, the system was solvated by adding SPC

ACCEPTED MANUSCRIPT (Single Point Charge) water. The starting structure was placed in a cubic box and it was centred with a distance of 0.9 Angstrom (Å) from the box boundary. The total charge of the system after solvation was noted; accordingly the system was neutralized by adding appropriate counter ions. Further, the solvated and neutralized system was energy minimized using steepest descent

T

algorithm for 50000 steps. However, the system converged well ahead within the specified steps,

IP

attaining a lowest potential energy. During these steps, the coulombs and Van der Waals

CR

interactions were maintained at a cut-off of 0.9 nm and 1.0 nm, respectively. Later, the system was equilibrated under NVT and NPT ensemble using Leap-frog integrator for 100 pico

US

seconds(ps) simultaneously maintaining Berendsen temperature coupling and Parinello-Rahman

AN

pressure coupling for 300k and 1 bar, respectively. The long range electrostatics was maintained by Particle-Mesh Ewald (PME) and short range electrostatics was maintained by coulombs and

M

Van der Waals at a cut off 1.0 nm each at a time-step of 10 femto seconds (fs) [36]. LINCS

ED

constraint algorithm [37] was applied to maintain the geometry of the molecules. Finally, a wellequilibrated system was then subjected to MD simulations for 30 nanoseconds (ns) each at 300 K

PT

without any constraints. The structural conformations were sampled at an interval of 2 ps. The

CE

final MD trajectories were analyzed for Root Mean Square Deviation (RMSD) of the protein backbone, Root Mean Square Fluctuations (RMSF) of the residues, Potential Energy distribution

AC

and radius of gyration (g_gyrate) during the production run. The trajectory plots for the above analysis were generated using XMGRACE software (http://plasma-gate.weizmann.ac.il/Grace/).

ACCEPTED MANUSCRIPT 2.7 Principal Component Analysis (PCA) and Free Energy Landscape (FEL) analysis for sampling the optimal conformations In order to address the collective motion of the macromolecules for both Apo [14] and Holo forms of RON4 with respect to its atomic coordinates and eigen vectors, principal

T

component analysis (PCA) was implemented in our study [38]. Since PCA is a multivariate

IP

statistical method, it linearly transforms and reduces the data to a positional covariance matrix

vectors. The covariance matrix is defined by

US

Ci=<(xi-)(xj -)>

CR

[39]. The 3N dimensional directions are derived by solving the covariance matrix of the eigen

AN

Where xi is a Cartesian coordinate of the ith Cα atom, N is the number of the Cα atoms considered, and represents the time average over all the configurations obtained in the

M

simulation. Here, for each data analysed, the first 10 Eigen vectors were generated, from which

ED

only first two principal component (PC) with minimum deviation and least cosine value of nearly 0.5 was considered for generating Gibb’s free energy using g_sham script. Further, free energy

PT

landscape analysis (FEL) was performed to sample the various conformations resulting from MD

CE

simulation. The PC1 and PC2 were considered as the two reaction coordinates, based on which

AC

2D representation of FEL was generated by implementing the equation: G(q1,q2)=-kB T ln P(q1,q2)

where KB is the Boltzmann constant, T is the simulation temperature, and P(q1,q2) is the normalized joint probability distribution [40]. Based on this, 2D projection of Free energy landscape was generated, and the minimum energy cluster for the macromolecule was determined thereby, suggesting the optimal conformations along the trajectory that can be utilized for further studies [41].

ACCEPTED MANUSCRIPT 2.8 Binding Free energy using Molecular Mechanics/Poisson--Boltzmann or generalized Born and surface area (MM/PBSA) The molecular mechanics energies combined with the Poisson-Boltzmann (MM/PBSA) method was used to calculate the binding free energies of small ligands to biological macromolecules.

T

GBINDING = Gcomplex- Gprotein-Gligand

IP

Where, Gcomplex represents the free energy of enzyme-substrate complex, and Gprotein and

CR

Gligand are the free energies of enzyme and substrate, respectively [42]. The solvation free

US

energy was calculated using an implicit-solvent model, in the MM/PBSA approach [43]. The solvation free energy is a combination of electrostatic terms, estimated by solving Poisson–

AN

Boltzmann (PB) equation and the non-polar solvation energy, separated as attractive and

M

repulsive interactions [44]. This was performed to decompose the energy determinants that plays

AC

CE

PT

ED

a key role in the establishment of intermolecular interactions during MD simulation.

ACCEPTED MANUSCRIPT 3. Results 3.1 RON4- Druggable binding site Prediction The druggable binding sites of RON4 were generated using SiteMap module of Schrodinger suite with a minimum of 15 site points per site. Among the binding sites generated

T

(Table 1 and Fig 1), site2 was identified as best druggable binding site based on SiteScore

IP

(1.046) and Druggability score (1.095). In addition to this, Site2 also had most significant scores

CR

in terms of its contact, hydrophobic, and hydrophilic components that corresponds to the plausible binding site (~1.0). Thus, based on this collective inference, site2 was finalized as

US

optimal druggable binding site for grid generation in virtual screening process. Moreover, this

AN

site was also found to be non-overlapping with AMA1 and  tubulin binding sites. The residues

M

that constitute the binding sites as per SiteMap are listed in Supplementary Table S1.

SiteScore Dscore volume

SITE1 SITE2 SITE3 SITE4 SITE5

1.026 1.046 0.997 0.993 0.988

2568.73 1400.47 1230.34 1080.11 608.14

AC

CE

1.056 1.095 1.019 1.017 1.001

exposure enclosure contact Hydro phobic 0.605 0.711 0.895 0.574 0.529 0.703 0.922 0.886 0.57 0.688 0.841 0.418 0.608 0.68 0.866 0.482 0.5 0.681 0.905 0.39

PT

Sites

ED

Table 1. Details of binding sites predicted by SiteMap Hydro philic 0.952 0.819 1.026 1.016 1.064

balance don/acc 0.603 1.081 0.407 0.474 0.367

0.841 1.191 1.794 1.392 1.074

ACCEPTED MANUSCRIPT 3.2 High Throughput Virtual Screening and identification of ophthalmic adaptive compounds based on PASS prediction analysis and MMGBSA Binding energy Out of 2,51,616 NCI compounds screened, only 2,29,865 compounds surpassed the initial Qikprop check. These compounds were further screened for drug likeness by Lipinski’s rule,

T

which resulted in 1,83,209 compounds. On subsequent filtering for excluding reactive functional

IP

groups, 1,46,670 compounds were identified and subjected to virtual screening against the

CR

potential druggable region (site2) of RON4. Among these 1,46,670 compounds, the top 10% compounds were scored as hits based on glide scoring and Hydrogen bonding analysis.

US

Subsequently, the resulting 14,667 compounds were fed to standard precision (SP) mode, which

AN

yielded 1,466 compounds (top 10% compounds of 14667). Finally, these 1,466 compounds were considered as best hits and were proceeded to highly stringent scoring by extra precision (XP),

M

which resulted in 147 compounds (top 10% of 1466 compounds). Among these 147 compounds,

ED

top 10 ranking compounds with higher Glide score (Table S2) were subjected to Prediction of Activity Spectra of Substances (PASS) analysis (Table S3). Initially, in the first tier of filtration

PT

for drug-likeness, Compound 3 (177021), Compound 5 (338690), Compound 9 (379447) among

CE

the top 10 hits lacked the drug-likeness property. Hence, the other seven compounds were considered for second tier of PASS prediction, which inferred that Compound 7 (286632) and

AC

Compound 10 (196524) to lack ophthalmic adaptive property. This overall filtration process resulted in 5 compounds that are ophthalmic adaptive: Compound 1 (338683), Compound 2 (405935), Compound 4 (159706), Compound 6 (686585), and Compound 8 (333739). Thus, these 5 compounds were re-ranked based on Molecular Mechanics Genralised Born Surface Area (MMGBSA) score (Table 2). It could be noted that Compound 1 (338683), Compound 8 (333739) to show significant binding free energy score of -91.728 Kcal/mol and -91.121

ACCEPTED MANUSCRIPT Kcal/mol, respectively, followed by other three compounds. Overall, the calculated binding free energy of the protein-ligand complexes ranged from −74 Kcal/mol to −91 Kcal/mol (Table S4). From this MMGBSA analysis, it could be inferred that Van der Waals (ΔGVdW) and solvation terms (ΔGsol) are the major contributors for favourable binding of ligand to RON4 protein. It

T

could also be noted that these two terms (ΔGVdW, ΔGsol) to act as significant factors in governing

IP

the binding affinity. To further prioritize, all these five compounds in complex with RON4 were

CR

subjected to MD simulation analysis in order to determine the stability of complex formation.

Druglikeness Property

Docking score

338683

Drug

-11.627

XP Gscore

AN

Compound 1

NSC Number

PASS Prediction

MMGBSA dG Bind (Kcal/mol)

-11.999

Corneal wound healing stimulator

-91.728

-10.292

-10.301

Anticataract, Antiglaucomic, Corneal wound healing stimulator, Ophthalmic drug

-91.121

M

Compound No.

US

Table 2. Top 5 ophthalmic adaptive compounds ranked based on MMGBSA Bind energy

333739

Drug

Compound 6

686585

Drug

-10.456

-10.456

Corneal wound healing stimulator

-86.258

Compound 4

159706

Drug

-11.116

-11.116

Corneal wound healing stimulator

-78.395

Compound 2

405935

Drug

-11.196

-11.222

Anticataract

-74.537

AC

CE

PT

ED

Compound 8

3.3 Density functional theoretical (DFT) analysis: The stability of interactions between protein-ligand complexes were determined based on the chemical reactivity of the small molecules based on theoretical energies using DFT. The theoretical energies HOMO, LUMO and gap energy of compounds are shown in Table 3. Further, the HOMO-LUMO plots were also generated for the compounds to analyse the atomic

ACCEPTED MANUSCRIPT contribution for these orbitals (Fig 2). In both the orbitals, red colour represents positive electron density and the blue colour denotes negative electron density. From the theoretical energies calculated, it could be noted that the gap energy for all the compounds were to be in average range of 0.5 Kcal/mol to 3.34 Kcal/mol, thereby signifying the molecular reactivity. The lowest

T

energy gap was observed in compound 4 (159706) with HLG of 0.5 Kcal/mol, as it had least

IP

electron donor ability (HOMO energy), yet had higher capability of accepting electrons (LUMO

CR

energy). Whereas, the highest energy gap was shown by compound 8 (333739) with 3.34 Kcal/mol, where the electron acceptor ability (LUMO energy) was lower when compared to its

US

electron donating (HOMO energy) ability. It is also noted that the HOMO energy is higher for all

AN

top five compounds when compared to the LUMO energy, thereby indicating a significant donor characteristic of the compounds leading to the strong stable interaction of compounds to RON4

M

protein.

PT

NSC Number 338683 333739 686585 159706 405935

HOMO energy (eV) -0.1864 -0.19305 -0.16524 -0.037364 -0.14223

LUMO energy (eV) -0.0764 -0.04809 -0.03197 -0.015235 -0.05597

HLG (Kcal/mol) 2.41 3.34 3.07 0.51 1.98

AC

CE

Compound No. Compound 1 Compound 8 Compound 6 Compound 4 Compound 2

ED

Table 3 Frontier orbital energies of top 5 compounds

3.4 Interaction analysis of the filtered Ophthalmic adaptive compounds The intermolecular interactions of the compounds shortlisted and validated using various filtering criteria such as PASS prediction, MMGBSA analysis and DFT analysis were tabulated (Table 4). From the results, it could be inferred that all the top five compounds were bound specifically to the site2 (active binding site) of RON4 as shown in Fig 3. In case of 338683, it

ACCEPTED MANUSCRIPT showed strong hydrogen bonds with Gln458, Asn481, Ala482, Gly488, Pro489, and Arg492 of RON4 (Fig 4a & Fig.S1a). Moreover, this complex was found to be stable, maintained by equal contribution of hydrogen bonded and hydrophobic interactions. In case of RON4-333739 complex, it formed three hydrogen bonds with Gln458, Thr491, Arg492 of RON4, it was also

T

stabilized by higher number of hydrophobic interactions (Fig 4b & Fig.S1b). RON4-686585

IP

complex was found to be stabilized by two hydrogen bonds (with Gly467, Arg492) and

CR

hydrophobic interactions (Fig 4c & Fig.S1c). In case of RON4-159706 complex, Asn481, Ala482, Gly488 residues of RON4 formed hydrogen bonds with the ligand stabilized by strong

US

hydrophobic interactions (Fig 4d & Fig.S1d). RON4-405935 complex showed (Fig 4e &

AN

Fig.S1e) three hydrogen bonds (formed by Gly488, Thr491, and Arg492 of RON4) and also showed some hydrophobic interactions. Among these top five compounds, the hydrophobic

M

interactions had higher impact on maintaining their stability, except for the RON4-338683

ED

complex, where the stability is contributed equally by hydrogen bonds. On overall analysis, it was noted that Arg492 to contribute for hydrogen bond in all the docked complexes except for

PT

RON4-159706. Apart from RON4-333739, all other docked complexes showed hydrophobic

AC

CE

interactions with Leu528.

ACCEPTED MANUSCRIPT Table 4. List of RON4 residues interacting with top 5 compounds

686585

Compound 4

159706

Compound 2

405935

T

Compound 6

IP

333739

CR

Compound 8

Gln458,Asn481, Trp473, Pro475, Leu490, Ala482,Gly488, Arg524, Ile525, Leu528 Pro489, Arg492 Asn481, Ala482, Ala486, Gln458,Thr491, Leu490, Val499, Val521, Arg492 Ile525 Leu463, Val468, Thr491, Pro498, Leu505, Val514, Gly467, Arg492 Leu528 Ala457,Leu463,Asn481, Asn481,Ala482, Ala482, Ala486, Leu490, Gly488 Leu528 Ile474, Asn481, Ala482, Gly488,Thr491, Ala486, Thr491, Pro498, Arg492 Leu528

US

338683

Hydrogen Bonds

M

Compound 1

Hydrophobic Interactions

AN

Compound No.

Residues involved in Interactions

NSC Number

ED

Furthermore, based on all the above validations, it could be inferred that all top five compounds shall act as the potential hits. However, to further prioritize, these docked complexes were

PT

subjected to essential molecular dynamics simulation studies, in order to validate the stability of

CE

complex formation at physiological conditions, and also to capture the conformational changes

AC

in the free energy landscape.

3.5 Molecular dynamics simulation analysis of docked complexes In order to compare the structural behaviour and stability of the top five compounds (based on MMGBSA and Ophthalmic adaptive), molecular dynamics studies was carried out. Further, MD trajectory analysis were performed and RMSD, RMSF, Radius of gyration, Hydrogen bond formation plots followed by PCA and FEL analysis was performed to assess the stability of complex formation at physiologically simulated conditions.

ACCEPTED MANUSCRIPT 3.5.1 RON4-338683 Complex From the RMSD plot of RON4-338683 Complex as shown in Fig S2a, it can be noted that the complex has maintained a maximum deviation of about ~1.0nm, after 150000ps not much deviation has been observed till the end of MD production run. On comparative analysis of

T

RMSD plots (Fig.S3b) of this Apo and Holo form of RON4, it is inferred Holo form has

IP

undergone maximum deviation than the Apo form with an average variation of about ~0.2nm,

CR

during 30ns of production run. Moreover, from the RMSF plot (Fig S2b), it could be inferred that residues Ser22 and Pro40 to have attained maximum flexibility of about ~2nm and ~3.8nm,

US

respectively. Other than these residues, the site2 (active binding site) residues has not showed

AN

any prominent flexibility. Fig S2c shows that, initially 7 hydrogen bonds were found to maintain the RON4-338683 complex, however on average only 4 hydrogen bonds has been observed to

M

constantly maintain the complex stability. The Radius of gyration plot as shown in Fig.S4a

ED

reveals that the compactness has been maintained at a range of ~3.5nm-3.3nm. The minimum potential energy structure of the complex (Fig.S5a) was observed at 17636th ps with potential

PT

energy of -2.9951e+06 kJ/mol, in this complex only Ala486, Asn454, Gln531, Arg524 were the

CE

RON4 residues that were found to be involved in the hydrogen bond formation with the compound (Fig.S5b).

AC

In addition, PCA analysis was carried out, thereby the overall combined motion of the Apo form and Holo form of RON4 was determined based on the projection of eigen vectors (Fig.S2d). The Holo form has split the frames and its projection represents that the Apo form has coordinated with the Holo form, and it is clearly observed that the Holo form had undergone conformation changes due to its interactions with the compound. Moreover, from the FEL analysis it has been inferred that the highly populated single minimum energy cluster and a least

ACCEPTED MANUSCRIPT populated cluster has been observed at 27120th ps (C1) and 9294th ps (C2), respectively (Fig.S2e). The minimum energy of the near native conformations of the docked complex at C1 was 0.0408 Kcal/mol. The representative structure of C1 cluster was extracted at 27120th ps and it was observed that complex has maintained it stability by one hydrogen bond formation

T

contributed by N481 (Fig.S2f & Fig.S7a, Table 5). Furthermore, the hydrophobic interactions

IP

made by Val493 and Salt Bridge by Arg524 were also found to enhance the stability of the

CR

docked complex. 3.5.2 RON4-333739 Complex

US

The RMSD plot of RON4-333739 Complex as shown in Fig 5a, reveals that the complex has

AN

increased deviation of ~1.2nm. Initially the complex has seemed to attain stability, yet it has not converged till the end of MD production run. On comparative analysis of the RMSD plots

M

(Fig.S3c) of Apo and Holo forms, it is observed that the Holo form has undergone minor

ED

deviation of less than ~0.1nm than the Apo form, yet both the forms have converged from 18000ps – 30000ps. The RMSF plot (Fig 5b) shows the residues have shown a maximum

PT

fluctuation of ~1.5nm. Fig 5c represents the intermolecular hydrogen bonds, in which 4-5

CE

hydrogen bonds were initially found to maintain the stability, yet on average only 4 hydrogen bonds were found to stabilize the complex. From the radius of gyration plot as shown in Fig.S4b,

AC

it could be inferred that the structural compactness to span in the range of ~3.6nm-3.3nm. The minimum potential energy structure of the complex (Fig.S5c) was extracted at 23194th ps with potential energy of -2.98783e+06 kJ/mol. This extracted complex was found to maintain its stability by hydrogen bonds formed by His487, Thr491 and Asn481 residues (Fig.S5d). Furthermore, PCA analysis resulted in overall combined motion of the Apo form and Holo form of RON4 based on the projection of eigen vectors (Fig 5d). The Holo form has completely

ACCEPTED MANUSCRIPT coordinated with its projection of Apo form, which reveals that the binding of compound to the RON4 has not altered its conformational space. In addition, the FEL analysis also substantiates the PCA analysis by providing a single highly populated minimum energy cluster obtained at 20442thps (C1) (Fig 5e). The minimum energy of this near native conformation of the docked

T

complex at C1 was 0.855 Kcal/mol. Moreover, from the representative structure of C1 cluster, it

IP

is noted that hydrogen bond formed by His487, Thr491 has maintained it stability. It is also

CR

inferred that hydrophobic interactions by RON4 residues Lys466, Trp473, Ala486, and Leu522 has the significant role in enhancing the stability of the docked complex (Fig 5f & Fig.S7b

US

Table 5). In addition, π-Stacking His487 with the compound was to found to promote the

M

3.5.3 RON4-686585 Complex

AN

binding affinity of the complex.

ED

The RMSD plot of RON4-686585 Complex (Fig.S8a), showed deviations at a range of ~0.8nm1.0nm. Initially, the complex has shown a constant increase in its deviation, however, it was

PT

found to stabilize over the last 10 ns (20000ps-30000ps). RMSD plots of the Apo and Holo forms

CE

of RON4, as shown in Fig.S3d revealed that the Holo form has undergone maximum deviation than the Apo form with a variation of about ~0.25nm. The residue Glu340 is the only residue

AC

with higher flexibility (~1.2nm), other than that all residues have fluctuated about ~0.25nm 0.95nm as shown in RMSF plot (Fig.S8b). The stability of the docked complex was found to be maintained by 6 hydrogen bonds, however, only 2-3 hydrogen bonds were observed across the production run (Fig.S8c). The Radius of gyration plot as shown in Fig.S4c reveals that the compactness has maintained a range of ~3.5nm-3.4nm. The minimum potential energy structure (-2.9954e+06 kJ/mol) was observed at 21818th ps, and its representative structure (Fig.S5e).

ACCEPTED MANUSCRIPT From the minimal energy complex, it could be inferred that Ala486, Ala457 to have majorly contributed to the hydrogen bonded interactions with the compound across the production run (Fig.S5f). In addition, from the PCA analysis, the projection of eigen vectors of RON4 both Apo and Holo

T

forms represents an overall combined motion (Fig.S8d). This plot reveals that the Holo form has

IP

not completely coordinated with its projection of Apo form, yet the motion has seen to be

CR

overlapped exactly on the compound binding site. This in turn reveals that the binding of this compound to RON4 has altered its conformational space drastically. Moreover, FEL analysis

US

provides a single highly populated minimum energy cluster (C1) (Fig.S8e) that corresponds to

AN

the energy value 0.82163 Kcal/mol. The representative structure of C1 cluster was extracted at 22224th ps and it was noted that the residues Gln458, Ala486, Glu495 has contributed to the

M

prominent hydrogen bond formation (Fig.S8f & Fig.S7c, Table 5). It could also inferred that

ED

hydrophobic interactions with Gln480, Pro489 of RON4 has contributed significantly for stable complex formation.

PT

3.5.4 RON4-159706 Complex

CE

From the RMSD plot of RON4-159706 Complex (Fig.S9a), protein backbone deviations at a range of ~0.8nm-1.0nm could be observed initially, however complex was found to stabilize

AC

over the last 20 ns (10000ps-30000ps). On Comparing the RMSD plots of the Apo and Holo forms, (Fig.S3e) it could be inferred that Holo form has undergone least deviation than the Apo form, nearly from 5000ps – 18000ps, the Holo form RMSD was found to be lower than that of Apo form. Later Holo form was found to maintain similar deviation pattern to that of Apo form and also found to get converged.

The RMSF plot (Fig.S9b) shows that the residues has

fluctuated about ~1.0nm, thereby, inferring that none of the residues has attained maximum

ACCEPTED MANUSCRIPT flexibility. The number of intermolecular hydrogen bonds in the docked complex was initially observed to be about 6-8, whereas, only 4 hydrogen bonds were found to be constantly maintained in the complex across the production run (Fig.S9c). The Radius of gyration plot as shown in Fig.S4d reveals that there was no prominent variation in the compactness of the

T

complex, as it has maintained ~3.5nm throughout the simulation of 30 ns. The minimum

IP

potential energy representative structure of the complex was extracted at 5796.5th ps with

CR

potential energy of -2.9943e+06 kJ/mol (Fig.S6a). This complex was found to be stabilized by hydrogen bonds formed by Ala504, Ala482, Asn481, and Gln531 to the compound (Fig.S6b).

US

Furthermore, in case of eigenvector projections, the Holo form was not completely found to be

AN

coordinated of Apo form, but the site 2 (active binding site) and its surrounding residues showed overlapping projections (Fig.S9d). This in turn reveals that the binding of this compound to the

M

RON4 has altered its conformational space. In addition, the FEL analysis also substantiates the

ED

PCA analysis by providing a varied conformation of the complex as 6 minimum energy clusters (Fig.S9e). Among the obtained clusters, the highly populated minimum energy cluster C1

PT

observed at 24776th ps was considered as the probable near native form. The representative

CE

structure of this docked complex (C1) was 0.3749 Kcal/mol. Moreover, from the representative structure of C1 cluster, it could be noted that the stability of

this docked complex was

AC

maintained by hydrogen bonds formed by Ile483, Pro489, Thr491, Ala504, Ala527 residues of RON4 (Fig.S9f & Fig.S7d, Table 5). In addition, the salt bridge formed by Arg535 and hydrophobic interactions of Ala456, Pro489 has significantly provided stability to the docked complex.

ACCEPTED MANUSCRIPT 3.5.5 RON4-405935 Complex The RMSD plot of RON4-405935 Complex as shown in Fig.S10a, infers that the complex has to show an initial deviation of ~0.9nm. However, the complex tend to attain the stability in the last 10ns (20000-30000ps) of MD production run. Comparative analysis on the RMSD plots

T

(Fig.S3f) of Apo and Holo forms showed no deviations. Moreover, both the forms were observed

IP

to show a similar pattern of convergence during the production run. The RMSF plot (Fig.S10b)

CR

shows the residues have attained a maximum fluctuation of ~1.0nm. The residues spanning the site2 (active binding site) has shown a slight flexibility. Fig.S10c which represents the

US

intermolecular hydrogen bonds of the docked complex reveals that 3-4 hydrogen bonds to be

AN

formed initially, yet an average only 2 hydrogen bonds has been observed to be constantly during the production run. From the radius of gyration plot as shown in Fig.S4e, it could be inferred that

M

the compactness has maintained a range of ~3.5nm-3.3nm. The minimum energy structure of the

ED

complex was extracted at 16200th ps with the potential energy of -2.9948e+06 kJ/mol (Fig.S6c). On analysis of this extracted complex, two hydrogen bonds formed by Gly488 and Ala486

PT

RON4 residues were found to confer stability to the complex (Fig.S6d).

CE

In addition, the projection of eigen vectors of Apo and Holo forms represents an overall combined motion. Fig.S10d reveals that the Holo form has not completely coordinated with the

AC

projection of Apo form, yet the motion has seen to be overlapped exactly on the compound binding site. This in turn reveals that the binding of this compound to alter the conformational space of RON4 in a significant manner. Moreover, the FEL analysis had shown two prominent minimum energy clusters (Fig.S10e). Among these clusters, the highly populated minimum energy cluster (C1) with 0.03415kcal/mol energy was considered for further analysis. The representative structure of C1 cluster was extracted at 21148th ps and it was noted that the

ACCEPTED MANUSCRIPT residues Gly488, Ile474 has contributed to the prominent hydrogen bond formation. It is also inferred that hydrophobic interactions by Trp473, Thr491 and Leu527 residues of RON4 to contribute to its complex stability (Fig.S10f & Fig.S7e, Table 5).

IP

T

Table 5. List of RON4 residues interacting with top 5 compounds as seen in the minimum energy cluster representative structures

Compound 6

686585

Compound 4

159706

Compound 2

405935

R524

K466,W473, A486, -H487, T491 L522 Q458, A486, -Q480, P489 E495 I483, P489, T491, R535 A456, P489 A504, 527

AN

333739

N481

M

Compound 8

V493, R524

ED

338683

Salt Bridges

Hydrogen Bonds

W473,T491, L527

G488, I474

--

π-Stacking -H487 ----

PT

Compound 1

Hydrophobic Interaction

US

NSC Number

CR

Residues involved in Interactions Compound No.

CE

3.6 MM/PBSA based Binding Free energy calculation and Decomposition analysis In this study, the binding free energy was also calculated based on MM/PBSA methodology,

AC

which provides a semi-quantitative estimation of inhibitor/ligand affinity to the protein. Table 6 represents the binding free energies and their components for the all the top 5 compounds. Among all the docked complexes, RON4-333739 complex was found to be most significant as it showed better binding affinity in terms of non-bonded terms Van der Waals energy (-250.570 ± 1.927 kJ/mol). The total binding free energy for the RON4-333739 complex was calculated as 179.122 ± 2.197 kJ/mol. While the total binding free energy of other complexes RON4-338683, RON4-159706 and RON4-405935 were calculated as -142.953 ± 1.728 kJ/mol , -147.930 ±

ACCEPTED MANUSCRIPT 1.885 kJ/mol and -145.017 ± 2.513 kJ/mol, respectively. Only RON4-686585 complex showed a least total binding free energy (-123.689 ± 2.743 kJ/mol) compared to other docked complexes. Furthermore, from the observations on the components of total binding free energy, it could be inferred that van der Waals and non-polar energies to be the major contributors that significantly

T

favours the binding of these compounds to RON4. In specific, van der Waals energy was found

IP

to dominate than the other components. Moreover, the binding free energies of all the

CR

compounds were negative values, thereby, indicating that these compounds to exhibit favourable

US

binding energetics.

ED

M

-38.005 ± 1.224

133.390 ± 2.286

-24.007 ± 0.149

-179.122 ± 2.197

-216.020 ± 3.173

-43.210 ± 2.878

156.736 ± 4.082

-21.369 ± 0.270

-123.689 ± 2.743

PT

-250.570 ± 1.927

CE

RON4338683 RON4333739 RON4686585 RON4159706 RON4405935

Van der Waal Electrostatic Polar solvation SASA energy Binding energy (kJ/mol) energy (kJ/mol) energy (kJ/mol) energy (kJ/mol) (kJ/mol) -239.035 ± 1.768 -62.473 ± 1.933 181.138 ± 2.367 -22.637 ± 0.158 -142.953 ± 1.728

-235.403 ± 2.488

-69.625 ± 2.161

179.493 ± 3.102

-22.348 ± 0.143

-147.930 ± 1.885

-206.157 ± 2.920

-0.912 ± 0.441

81.807 ± 2.134

-19.698 ± 0.165

-145.017 ± 2.513

AC

Complex

AN

Table 6 Decomposed Components of the binding free energy (kJ/mol) for the RON4 complexes with top five compounds as determined during MM/PBSA calculation

Per-residue free energy contributions of complexes In order to characterize and identify the key residues of the compounds and RON4 in all the docked complexes, per-residue free energy decompositions was performed to obtain their individual residue energy contributions, as shown in Fig6. On overall analysis of binding free energies of all the residues, it was inferred that the residue Asp399 favoured the binding of

ACCEPTED MANUSCRIPT compound 686585 to RON4, while in other docked complexes it showed unfavourable interactions. Furthermore, it was observed that among the site2 (active binding site), the residues Leu490, Pro494, Leu505, Ile525, Leu528, Ala457, Val468, Trp473, Lys474, Met475, Ala482, Lys485, Ala486, Pro489, Val493 were observed to play a crucial role in the binding of all the

T

compounds to RON4. On cumulative analysis of the MMPBSA decomposition, binding free

CR

IP

energy of RON4-333739 complex was found to be highly significant when compared to others.

4. Conclusion

US

In this study, the compounds 338683, 333739, 686585, 159706, 405935 has been predicted as a

AN

potent ophthalmic adaptive inhibitors for targeting RON4 of T.gondii through virtual screening, PASS prediction and MMGBSA analysis. Moreover, the DFT studies revealed that 159706,

M

405935 to be highly reactive due to its lowest HOMO-LUMO energy gap compared to

ED

compounds. However, among the identified compounds, the compounds 333739, 159706, 405935 were shortlisted as highly potential based on molecular docking, molecular dynamics

PT

simulation studies, as these compounds formed significant non-bonded interactions with the

CE

druggable cavity of RON4, thereby, forming stable protein-drug complex. Furthermore, from the projection of eigenvectors, it could be noted that 333739 to be the only compound to confer

AC

significant arrest of conformational shift in RON4. In addition, the MM/PBSA analysis revealed that the compound 333739 to have lowest binding free energy, which was significantly favoured by Van der Waals terms. Thus, by this study we propose 333739 (benzyl 2-((2-(((benzyloxy) carbonyl) amino)-3-hydroxybutanoyl) amino) propanoate) to be the most potential inhibitor Toxoplasma gondii RON4, which shall prove to be effective in the treatment of ocular

ACCEPTED MANUSCRIPT toxoplasmosis. However, further experimental studies are required to validate the efficacy of this molecule.

Acknowledgement

T

This study was supported by the Department of Biotechnology (DBT), Government of India,

IP

under Rapid Grant for Young Investigator (RGYI) scheme [BT/PR6476/GBD/27/496/2013,

AC

CE

PT

ED

M

AN

US

CR

05/09/2013].

ACCEPTED MANUSCRIPT References

AC

CE

PT

ED

M

AN

US

CR

IP

T

[1] J.G. Montoya, Laboratory diagnosis of Toxoplasma gondii infection and toxoplasmosis, The Journal of infectious diseases 185 Suppl 1 (2002) S73-82. [2] Indhuja Thirumudi, Umashankar Vetrivel, Mahalakshmi B, Lily Theresa K, MadhavanHN, Insights on drug targeting of Toxoplasma gondii host invasion protein: A review 8 (3) (2015). [3] J.A. Kovacs, Efficacy of atovaquone in treatment of toxoplasmosis in patients with AIDS. The NIAIDClinical Center Intramural AIDS Program, Lancet (London, England) 340 (8820) (1992) 637–638. [4] D.C. McFadden, S. Tomavo, E.A. Berry, J.C. Boothroyd, Characterization of cytochrome b from Toxoplasma gondii and Q(o) domain mutations as a mechanism of atovaquone-resistance, Molecular and biochemical parasitology 108 (1) (2000) 1–12. [5] J.C. Boothroyd, Toxoplasma gondii: 25 years and 25 major advances for the field, International journal for parasitology 39 (8) (2009) 935–946. [6] S. Besteiro, J.-F. Dubremetz, M. Lebrun, The moving junction of apicomplexan parasites: a key structure for invasion, Cellular microbiology 13 (6) (2011) 797–805. [7] L.M. Weiss, A. Fiser, R.H. Angeletti, K. Kim, Toxoplasma gondii proteomics, Expert review of proteomics 6 (3) (2009) 303–313. [8] P.J. Bradley, C. Ward, S.J. Cheng, D.L. Alexander, S. Coller, G.H. Coombs, J.D. Dunn, D.J. Ferguson, S.J. Sanderson, J.M. Wastling, J.C. Boothroyd, Proteomic analysis of rhoptry organelles reveals many novel constituents for host-parasite interactions in Toxoplasma gondii, The Journal of biological chemistry 280 (40) (2005) 34245–34258. [9] S. Besteiro, A. Michelin, J. Poncet, J.-F. Dubremetz, M. Lebrun, Export of a Toxoplasma gondii rhoptry neck protein complex at the host cell membrane to form the moving junction during invasion, PLoS pathogens 5 (2) (2009) e1000309. [10] H. Takemae, K. Kobayashi, T. Sugi, Y. Han, H. Gong, A. Ishiwa, F.C. Recuenco, F. Murakoshi, R. Takano, Y. Murata, K. Nagamune, T. Horimoto, H. Akashi, K. Kato, Toxoplasma gondii RON4 binds to heparan sulfate on the host cell surface, Parasitology international 67 (2) (2018) 123–130. [11] H. Takemae, T. Sugi, K. Kobayashi, H. Gong, A. Ishiwa, F.C. Recuenco, F. Murakoshi, T. Iwanaga, A. Inomata, T. Horimoto, H. Akashi, K. Kato, Characterization of the interaction between Toxoplasma gondii rhoptry neck protein 4 and host cellular β-tubulin, Scientific reports 3 (2013) 3199. [12] D.L. Alexander, S. Arastu-Kapur, J.-F. Dubremetz, J.C. Boothroyd, Plasmodium falciparum AMA1 binds a rhoptry neck protein homologous to TgRON4, a component of the moving junction in Toxoplasma gondii, Eukaryotic cell 5 (7) (2006) 1169–1173. [13] M. Lebrun, A. Michelin, H. El Hajj, J. Poncet, P.J. Bradley, H. Vial, J.F. Dubremetz, The rhoptry neck protein RON4 re-localizes at the moving junction during Toxoplasma gondii invasion, Cellular microbiology 7 (12) (2005) 1823–1833. [14] U. Vetrivel, H. Nagarajan, I. Thirumudi, Design of inhibitory peptide targeting Toxoplasma gondii RON4-human β-tubulin interactions by implementing structural bioinformatics methods, Journal of cellular biochemistry 119 (4) (2018) 3236–3246. [15] S.G. V.Umashankar, Drug Discovery: An appraisal, International Journal of Pharmacy and Pharmaceutical Sciences (4) (2015).

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

[16] S.G. Vetrivel Umashankar, Chemoinformatics and its Applications, General, Applied and Systems Toxicology (2011). [17] Z. Guo, B. Li, L.-T. Cheng, S. Zhou, J.A. McCammon, J. Che, Identification of protein-ligand binding sites by the level-set variational implicit-solvent approach, Journal of chemical theory and computation 11 (2) (2015) 753–765. [18] T.A. Halgren, Identifying and characterizing binding sites and assessing druggability, Journal of chemical information and modeling 49 (2) (2009) 377–389. [19] G.M. Sastry, M. Adzhigirey, T. Day, R. Annabhimoju, W. Sherman, Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments, Journal of computer-aided molecular design 27 (3) (2013) 221–234. [20] J.H. Voigt, B. Bienfait, S. Wang, M.C. Nicklaus, Comparison of the NCI open database with seven large chemical structural databases, Journal of chemical information and computer sciences 41 (3) (2001) 702–712. [21] W.-D. Ihlenfeldt, J.H. Voigt, B. Bienfait, F. Oellien, M.C. Nicklaus, Enhanced CACTVS browser of the Open NCI Database, Journal of chemical information and computer sciences 42 (1) (2002) 46–57. [22] A. John, V. Umashankar, A. Samdani, M. Sangeetha, S. Krishnakumar, P.R. Deepa, In Silico Structure Prediction of Human Fatty Acid Synthase-Dehydratase: A Plausible Model for Understanding Active Site Interactions, Bioinformatics and biology insights 10 (2016) 143–154. [23] A. John, M. Sivashanmugam, V. Umashankar, S.K. Natarajan, Virtual screening, molecular dynamics, and binding free energy calculations on human carbonic anhydrase IX catalytic domain for deciphering potential leads, Journal of biomolecular structure & dynamics 35 (10) (2017) 2155– 2168. [24] A. John, U. Vetrivel, K. Subramanian, P.R. Deepa, Comparative docking of dual conformations in human fatty acid synthase thioesterase domain reveals potential binding cavity for virtual screening of ligands, Journal of biomolecular structure & dynamics 35 (6) (2017) 1350–1366. [25] M. Sivashanmugam, C. Raghunath, U. Vetrivel, Virtual screening studies reveal linarin as a potential natural inhibitor targeting CDK4 in retinoblastoma, Journal of pharmacology & pharmacotherapeutics 4 (4) (2013) 256–264. [26] S. Parasuraman, Prediction of activity spectra for substances, Journal of pharmacology & pharmacotherapeutics 2 (1) (2011) 52–53. [27] V.V. Poroikov, D.A. Filimonov, W.-D. Ihlenfeldt, T.A. Gloriozova, A.A. Lagunin, Y.V. Borodina, A.V. Stepanchikova, M.C. Nicklaus, PASS biological activity spectrum predictions in the enhanced open NCI database browser, Journal of chemical information and computer sciences 43 (1) (2003) 228– 236. [28] U. Vetrivel, S. Muralikumar, B. Mahalakshmi, K. Lily Therese, H.N. Madhavan, M. Alameen, I. Thirumudi, Multilevel Precision-Based Rational Design of Chemical Inhibitors Targeting the Hydrophobic Cleft of Toxoplasma gondii Apical Membrane Antigen 1 (AMA1), Genomics & informatics 14 (2) (2016) 53–61. [29] S. Muralikumar, U. Vetrivel, A. Narayanasamy, U. N Das, Probing the intermolecular interactions of PPARγ-LBD with polyunsaturated fatty acids and their anti-inflammatory metabolites to infer most potential binding moieties, Lipids in health and disease 16 (1) (2017) 17.

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

[30] T. Hou, J. Wang, Y. Li, W. Wang, Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations, Journal of chemical information and modeling 51 (1) (2011) 69–82. [31] S. Salentin, S. Schreiber, V.J. Haupt, M.F. Adasme, M. Schroeder, PLIP: fully automated proteinligand interaction profiler, Nucleic acids research 43 (W1) (2015) W443-7. [32] A.D. Bochevarov, E. Harder, T.F. Hughes, J.R. Greenwood, D.A. Braden, D.M. Philipp, D. Rinaldo, M.D. Halls, J. Zhang, R.A. Friesner, Jaguar: A high-performance quantum chemistry software program with strengths in life and materials sciences, Int. J. Quantum Chem. 113 (18) (2013) 2110– 2142. [33] M. Sivashanmugam, S. K N, U. V, Virtual screening of natural inhibitors targeting ornithine decarboxylase with pharmacophore scaffolding of DFMO and validation by molecular dynamics simulation studies, Journal of biomolecular structure & dynamics (2018) 1–15. [34] M.J. Abraham, T. Murtola, R. Schulz, S. Páll, J.C. Smith, B. Hess, E. Lindahl, GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX 1-2 (2015) 19–25. [35] D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A.E. Mark, H.J.C. Berendsen, GROMACS: fast, flexible, and free, Journal of computational chemistry 26 (16) (2005) 1701–1718. [36] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, Molecular dynamics with coupling to an external bath, The Journal of Chemical Physics 81 (8) (1984) 3684–3690. [37] B. Hess, H. Bekker, H.J.C. Berendsen, Fraaije, Johannes G. E. M., LINCS: A linear constraint solver for molecular simulations, J. Comput. Chem. 18 (12) (1997) 1463–1472. [38] C.C. David, D.J. Jacobs, Principal component analysis: a method for determining the essential dynamics of proteins, Methods in molecular biology (Clifton, N.J.) 1084 (2014) 193–226. [39] A. Amadei, A.B. Linssen, H.J. Berendsen, Essential dynamics of proteins, Proteins 17 (4) (1993) 412– 425. [40] E. Papaleo, P. Mereghetti, P. Fantucci, R. Grandori, L. de Gioia, Free-energy landscape, principal component analysis, and structural clustering to identify representative conformations from molecular dynamics simulations: the myoglobin case, Journal of molecular graphics & modelling 27 (8) (2009) 889–899. [41] H. Nagarajan, U. Vetrivel, Demystifying the pH dependent conformational changes of human heparanase pertaining to structure-function relationships: an in silico approach, Journal of computer-aided molecular design (2018). [42] R. Kumari, R. Kumar, A. Lynn, g_mmpbsa--a GROMACS tool for high-throughput MM-PBSA calculations, Journal of chemical information and modeling 54 (7) (2014) 1951–1962. [43] R.C.C. M, H.-P. Hsieh, M.S. Coumar, Delineating the active site architecture of G9a lysine methyltransferase through substrate and inhibitor binding mode analysis: a molecular dynamics study, Journal of biomolecular structure & dynamics (2018) 1–36. [44] S. Genheden, U. Ryde, The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities, Expert opinion on drug discovery 10 (5) (2015) 449–461.

ACCEPTED MANUSCRIPT Figure Legends

IP

T

Fig 1 Surface representation of the RON4 protein along with the predicted binding sites obtained from SiteMap predictions The colour representation RON4 (teal), Site1 (red), Site2 (yellow), Site3 (orange), Site4 (purple), Site5 (pale pink)

US

CR

Fig 2 HOMO-LUMO plots of (a) 338683 HOMO-LUMO (b) 333739 HOMO-LUMO (c) 686585 HOMO-LUMO (d) 159706 HOMO-LUMO (e) 405935 HOMO-LUMO; the red colour represents positive electron density and the blue colour represents negative electron density.

AN

Fig 3 Surface representation of the top5 compounds docked with RON4 (a) Represents RON4 (teal) with all ligands (stick representation) docked at active site (site 2 –yellow); (b) 338683-Red (c) 333739-Blue (d) 686585-Green (e) 159706-Magenta (f) 405935- Black

ED

M

Fig 4 2D interaction plots of top 5 compounds shortlisted based on the MM-GBSA score and ophthalmic adaptive activity with RON4 (a) 338683 (b) 333739 (c) 686585 (d) 159706 (e) 405935

AC

CE

PT

Fig 5 MD simulation analysis of RON4-333739 Complex (a) RMSD plot showing the deviation of ~1.2nm (b) RMSF plot representing the fluctuations of RON4 residues in the docked complex (c) Number of hydrogen bonds observed during the production run ; PCA and FEL analysis of RON4-333739 Complex (d) Projection of Eigen vectors of both Apo (orange) and Holo (blue) form of RON4 (e) 2D projection of the Free energy landscape contour map (f) Interaction of minimum energy cluster representative structure of RON4333739 complex observed 20442thps Fig 6 Per-residue free energy contributions of complexes (a) RON4-338683 complex (b) RON4-333739 complex (c) RON4-686585 complex (d) RON4-159706 complex (e) RON4405935 complex

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6