A Computational Approach to the Study of the Binding Mode of S1P1R

structural analogues synthesized in-house and 1145 known S1P1R agonists. ...... Jo, S.; Klauda, J. B.; Im, W.; Lim, J. B. CHARMM-GUI Membrane Buil...
0 downloads 0 Views 1MB Size
Subscriber access provided by UNIV OF TEXAS DALLAS

Pharmaceutical Modeling

A Computational Approach to the Study of the Binding Mode of S1P1R Agonists Based on the Active-like Receptor Model Yonghui Chen, Tianqi Liu, Qiumu Xi, Wenqiang Jia, Dali Yin, and Xiaojian Wang J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00764 • Publication Date (Web): 11 Mar 2019 Downloaded from http://pubs.acs.org on March 12, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

A Computational Approach to the Study of the Binding Mode of S1P1R Agonists Based on the Active-like Receptor Model Yonghui Chena, Tianqi Liub, Qiumu Xib, Wenqiang Jiaa, Dali Yina,b, Xiaojian Wang*a,b aState

Key Laboratory of Bioactive Substances and Functions of Natural Medicines,

Institute of Materia Medica, Peking Union Medical College and Chinese Academy of Medical Sciences, Beijing 100050, P.R. China. bDepartment

of Medicinal Chemistry, Beijing Key Laboratory of Active Substances

Discovery and Druggability Evaluation, Institute of Materia Medica, Peking Union Medical College and Chinese Academy of Medical Sciences, Beijing 100050, P.R. China.

ABSTRACT: Sphingosine-1-phosphate receptor 1 (S1P1R), a member of the G protein-coupled receptor (GPCR) family, is an attractive protein target for the treatment of autoimmune diseases, and a diverse array of S1P1R agonists have been developed. Rational drug design based on S1P1R remains challenging due to the limited information available on the binding mode between S1P1R and its agonists. In this work, the active-like state of S1P1R was modeled via Gaussian accelerated molecular dynamics (GaMD) based on its inactive form, which was further validated by docking 1 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 35

studies with two representative S1P1R agonists. Moreover, with the usage of the induced active-like state, the binding mode between S1P1R and its agonists was studied through molecular dynamics simulations and MM-GBSA calculations. The results of those studies indicated that four groups of binding site residues were the major contributors to the ligand and receptor interactions. In addition, this model was verified by five chemically similar compounds synthesized in-house and 1145 known S1P1R agonists collected from the BindingDB database. The elucidation of the key binding characteristics will further complete the cognition of S1P1R, which can guide the rational design of novel S1P1R agonists.

INTRODUCTION Spingosine-1-phosphate receptor 1 (S1P1R), a G protein-coupled receptor (GPCR), has emerged as an attractive drug target.1-3 The activation of S1P1R by agonists leads to sustained internalization and degradation of the receptor, which further results in immunosuppression, providing a potential platform for treating autoimmune diseases such as plaque psoriasis, lupus, relapsing-remitting multiple sclerosis (RRMS).4, 5 Known agonists of S1P1R (Figure 1) can be classified into three basic types: class I ligands are lipid-like compounds that mimic the interaction between S1P1R and the sphingophospholipid headgroup of S1P, as exemplified by S-FTY720-P (active phosphorylated fingolimod, 1);6 class II agonists, such as ASP4058 (2), do not require polar headgroups and instead form specific binding interactions through their hydrophobic tail;7,

8

class III agonists, APD334 (3) for example, join the polar

headgroup and the hydrophobic tail together with a linker.9 2 ACS Paragon Plus Environment

Page 3 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

HO O HO P O

5

CH3

N

HO

O OH

CF3

S-FTY720-P 1

CF3

HN O

O

N

N H

NH2

CF3

H

N O

APD334 3

ASP4058 2

Cl

O HO

N N

BAF312 4

O

N

HO CF3

O

N

N O

N

GSK2018682 5

N

H N

O

N

HO N

O

O

RPC1063 6

Figure 1. Structures of representative S1P1R agonists. Notably, despite the structural diversity of S1P1R agonists, a proximal disubstituted aromatic ring is frequently observed in the hydrophobic tail of these molecules (Figure 1).10 We envisioned that this moiety forms particular interactions with corresponding residues of S1P1R. Therefore, it is essential to identify the binding mode between S1P1R and its agonists, and this would provide sufficient guidance for the further development of S1P1R agonists. However, only the crystal structure of S1P1R bound to its selective antagonist sphingolipid mimic (PDB code: 3V2Y and 3V2W) has been reported, and this structure does not provide sufficient information for rationally assessing the binding characteristics of S1P1R with its agonists.11 In this work, we reported the active-like state of the S1P1R model based on its inactive state using GaMD.12 Comparing to the inactive state, the active-like model was observed with the outward displacement of the TM5 and TM6 in the cytoplasmic end, which resulted in the breaking of the “ion lock” between the Arg78 and Glu250. Meanwhile, the residue Trp269 as a “toggle-switch” was rotated almost 90 degree. Subsequent docking studies were conducted to determine the accuracy of the model. 3 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 35

Furthermore, a series of conventional molecular dynamics (cMD) simulations was performed on the active-like S1P1R in complex with three clinical drug candidates, BAF312 (4), GSK2018682 (5), and RPC1063 (6).13-15 The effective binding energies were decomposed into the contributions from individual residues using the MM-GBSA energy decomposition scheme, which led to the observation of four groups of binding site residues in S1P1R.16, 17 In addition, computational and experimental approaches were used to examine the involvement of the binding site residues in the interaction between S1P1R and its agonists.

MATERIALS AND METHODS System preparation. The crystal structure of the inactive state of human S1P1R (PDB code: 3V2Y) was used for molecular docking and simulations. The antagonist ML056 and the T4lysozyme were removed, and the missing residues and intracellular loops were repaired using the protein preparing tools in Discovery Studio 4.5. In particular, the sequences of the missing loop (ICL3) were obtained from the UniProt database (UniProt ID: P21453).18 After that, the S1P1R agonist in clinical trials, compound 3, with a carboxylic acid as the polar headgroup and a disubstituted benzene as the hydrophobic tail, was docked into the receptor using the Glide program (Schrödinger 2009)19 with the standard precision (SP) setting. During the docking work, the Protein Preparation Wizard was used to add hydrogen atoms; fix bond orders; assign partial charges with the OPLS force field. Meanwhile, the LigPrep component was applied to pretreat the 4 ACS Paragon Plus Environment

Page 5 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

ligand with correction of the structure, conformation generation, especially the ionization of the polar head of compound 3. After that, an autosized binding box of the protein centered on the antagonist position was generated. The van der Waals radius scaling factor were set to 1.0 and the partial atomic charge cutoff were 0.25. Finally, the rational docking pose with the top Gscore was selected for subsequent simulations. Before the dynamic simulation, the CHARMM-GUI server was used to embed the receptor-agonist complexes in a water-lipid box with dimensions of 79.79 ×79.83 × 109.13

3

and neutralized with potassium and chloride ions (Figure 2).20 The lipid

bilayer was composed of 80% 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) and 20% cholesterol.

Figure 2. Overview of the initial inactive S1P1R/lipid bilayer/compound 3/water/KCl system  Conventional and Gaussian accelerated molecular dynamics.

5 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 35

On the basis of the initial system described above, the molecular dynamics was evaluated by means of Amber14 software with the leaprc.ff14SB and leaprc.lipid14 force field parameters.21, 22 A 20 ns cMD simulation was performed in the isothermalisobaric (NPT) ensemble at a simulation temperature of 303 K to ensure the system reached equilibrium. The Verlet integrator with an integration time step of 0.002 ps and SHAKE constraints for covalent bonds involving hydrogen atoms was applied to the simulations. For electrostatic interactions, atom-based truncation was carried out using the particle mesh Ewald method, and the switch van der Waals function with a 2.00 nm cutoff was used for atom-pair lists. The final frame of the trajectory was extracted and subjected to three independent GaMD. The GaMD is a biomolecular enhanced sampling method implemented by adding a harmonic boost potential to smoothen the system potential energy surface and has been implemented in Amber 14 package, which was reported by Miao et al. in previous study. This method has been applied to alanine dipeptide, chignolin, lysozyme, HIV protease. 23-25

In our study, the GaMD simulations were used to simulate the conformational dynamics of S1P1R. Bonds containing hydrogen atoms were restrained with the SHAKE algorithm, and the bond interactions involving H-atoms were omitted. Relative geometrical tolerance for coordinate resetting shake was set by 0.000001 Å, and the nonbonded cutoff was 8.00 Å. The Langevin thermostat was also used to maintain a temperature of 303 K and the simulation was subjected with isothermal constantvolume. A 2 ns short cMD simulation was processed to collect the maximum, minimum, 6 ACS Paragon Plus Environment

Page 7 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

average, and standard deviation values of the system potential for calculating the GaMD acceleration parameters, and a 6 ns biasing molecular dynamics simulation was performed after adding the boost potential. Finally, the GaMD production was carried with dual boost on both dihedral and total potential energy, and the threshold energy was set to the lower bound (E=Vmax). The average and standard deviation of the system potential energies were calculated every 100 ps. The upper limit of the boost potential standard deviation of the dihedral and total potential energies were both set to 6.0 kcal/mol. Three independent GaMD simulations (One for 900 ns and two for 600 ns) were performed with random atomic velocity initializations. Principal component analysis (PCA) and reweighted free energy calculations. After the GaMD simulations, the calculation of the principal component analysis (PCA) describing the direction and amplitude of the receptor motions was performed with AmberTool17 employing the Cpptraj component.26 The motion of the amino acids Ile259, Ile260, and Val261 on TM6 of S1P1R was analyzed since these amino acids were related to the outward displacement of TM6 in the cytoplasmic end, and this displacement was one of the microscopic characteristics of GPCR activation. The top two principal components, PC1 and PC2, representing the largest contributors to the overall fluctuation of these three amino acids (50−60%), provided the reaction coordinates for potential mean force (PMF) calculations. The PyReweighting toolkit was used to reweight the GaMD simulations to calculate the two-dimensional (2D) PMF profile and to examine the boost potential distribution with a unit bin size of 0.01 Å.27 Moreover, the RMSD values for the simulation trajectory relative to the inactive 7 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 35

state of the S1P1R structure were also calculated to obtain the one-dimensional (1D) PMF profile with a unit size of 0.1 Å. The snapshot of the active-like state of S1P1R could be obtained by searching in the low-energy sites of the two PMF profiles. Molecular docking with representative S1P1R agonists. Since a potential active-like state of S1P1R was extracted from the GaMD simulation, two S1P1R agonists, compounds 1 and 2, were docked into the active site of S1P1R using the Glide component implemented in Schrodinger 2009 to verify the rationality of the model. For the protein structure, the Protein Preparation Wizard was used to remove lipid bilayers, water molecules, and ions; add hydrogen atoms; fix bond orders; assign partial charges with the OPLS force field; and the structures of the two ligands were processed by the LigPrep component. After this, an autosized binding box of the protein centered on docked agonist compound 3 was generated using the Receptor Grid Generation component of Glide. The van der Waals radius scaling factor and the partial atomic charge cutoff were set to 1.0 and 0.25, respectively. All the prepared structures of the two ligands were docked into the active-like state of S1P1R and scored by the SP scoring mode. Molecular dynamics simulations and MM-GBSA calculations. In a further study, the MM-GBSA method was used to investigate the energetic contributions of each residue surrounding the binding pocket. The basic concept of the MM-GBSA approach is that the free energy of binding can be obtained by calculating the end points of the thermodynamic cycle of ligand binding. The binding free energy, ΔGbind, is calculated as follows: 8 ACS Paragon Plus Environment

Page 9 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

ΔGbind = ΔGcomplex - ΔGprotein - ΔGligand

(1)

ΔGbind = ΔH-TΔS

(2)

Each term is estimated as equation (7): ΔGbind = ΔGMM + ΔGsolv - TΔS

(3)

where ΔGMM is the molecular mechanics free energy, ΔGsolv is the solvation free energy, and TΔS is the entropy term. The ΔGMM can be further divided into the van der Waals contribution (ΔEvdw) and the electrostatic interaction energies (ΔEele). In addition, the solvation free energy, ΔGsolv, is composed of polar and nonpolar components: ΔGMM = ΔEvdw + ΔEele

(4)

ΔGsolv = ΔGnonpol solv +ΔGele solv

(5)

The contribution of the polar solvation energy, ΔGele

solv,

is calculated with the

generalized Born (GB) implicit solvent model, whereas the nonpolar part of the solvation energy, ΔGnonpol solv, is dependent on the solvent-accessible surface area (SA). TΔS is the contribution arising from changes in the degrees of freedom of the solute molecules, and this is not considered here. Thus, the values reported for the MM-GBSA calculations should be considered “effective energies” rather than free energies.28 Initially, three S1P1R agonists, compounds 4, 5, and 6, were docked into the active site of S1P1R. The subsequent 40 ns conventional molecular dynamics simulations were performed with the same steps as the cMD mentioned above. The time-dependence of the RMSD of the receptor backbone atoms was calculated for each trajectory using the Cpptraj component of AmberTool17. 9 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 35

The decomposed binding free energy was calculated for each 40 ns cMD trajectory to determine the binding site residues through the MM-GBSA energy decomposition method implemented in Amber14. Snapshots at 10 ps intervals were extracted from the last 10 ns of the cMD trajectories, and the decomposed binding free energies were averaged over the ensemble of conformers produced. The binding interaction energies of the individual residues included four terms: the van der Waals contribution (ΔEvdw), the electrostatic contribution (ΔEele), and the two solvation contribution terms (ΔGpolar solv +

ΔGnon-polar solv).

Synthesis and biological evaluation of oxadiazole derivatives. Three compounds (14a-14c, Figure 3) were selected from Liu et al’s manuscript,29 and the compounds with modification on polar head (14d) and hydrophobic tail (13e) were synthesized and subjected to biological evaluations to examine the predicted binding mode of the active-like S1P1R with its agonists. The synthetic pathway utilized for the preparation of the oxadiazole derivatives was outlined in Scheme 1. Among these compounds, 14a and 13e were used to reveal the interactions between the ligand’s polar headgroup and the hydrophilic subpocket of the receptor, while 14a-d were applied to elucidate the interactions between the terminal disubstituted benzene of the ligand and the hydrophobic subpocket of the receptor. These compounds were also docked into the active-like state of S1P1R using the Glide component and scored by the SP scoring mode with default sets.

10 ACS Paragon Plus Environment

Page 11 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

R

O N

O

14a 14b 14c

Cl

N

H N

R=CF3 R=CH3 R=H

COOH

Figure 3. The structures of the compound 14a, 14b and 14c. Cl

Cl

a

NC

Br

NC

7 F3 C

O

Cl

c

HO N

O

NC

CHO

O

8 COOH

R1

Cl

b

O

H2 N

9

F3 C

O N

R1

F3 C

e

Cl

10

N

R1

d

O N

Cl

N

O

CHO

O 11a 11d

f

F3 C R1

R1=isopropoxy R1=methoxyl

O N N

13d 13e

12a 12d

g,h

Cl H N

F3 C R1

R2 '

R1=methoxyl, R2'=CH2COOCH3 R1=isopropoxy, R2'=CH2CH2CH3

14d

R1=isopropoxy R1=methoxyl

O N N

Cl H N

R2

R1=methoxyl, R2=CH2COOH

Scheme 1. Synthesis of oxadiazole derivatives. Reagents and conditions: (a) i-PrMgCl, N-formylpiperidine, THF, 0°C, 4 h; (b) Ethylene glycol, PTSA, TOL, reflux, 18 h; (c) Hydroxylamine hydrochloride, NaHCO3, MeOH, reflux, 5 h; (d) HOBt, EDCI, K2CO3, DMF, 110°C, 7 h; (e) 2N HCl, acetone, 60°C, 5 h; (f) Amino acid ester or propylamine, DIPEA, AcOH, NaCNBH3, MeOH, CH2Cl2, r.t., 5 h; (g) 1 N LiOH, MeOH, reflux, 5 h. (h) 0.5 N HCl, acidification.

RESULTS AND DISCUSSION Activation of S1P1R observed in GaMD simulation. In the previous study, the GaMD method enabled direct observation of the activation of the GPCR and identification of the active-like and inactive states.30 It was well acknowledged that the activation of GPCR was usually characterized by three features: (i) a toggle-switch of the Trp residue on TM6;31, 32 (ii) the breaking of an “ionic lock” 11 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 35

in the cytoplasmic end;33 (iii) and the outward displacement of TM5 and TM6 in the cytoplasmic end by 6~7 Å.34, 35 These features drove the inactive receptor into an active state that was similar to the X-ray structure of opsin.36, 37 In our simulation, the potential of mean force (PMF) maps of the three independent trajectories were constructed by reweighting the GaMD simulation according to the applied boost potential (Figure S1). The 900 ns GaMD simulation was chosen for further study because the low-energy conformations of the receptor were then observed clearly in the stable state, which showed the inactive and active-like conformations of S1P1R. The one-dimensional (1D) PMF was calculated for the RMSD values of the S1P1R backbone atoms related to the initial structure (Figure 4A), suggesting that the system was situated in a low-energy state with an RMSD value greater than 4 Å. In addition, the two-dimensional (2D) PMF was calculated for the principal components PC1 vs. PC2 (Figure 4B). Then, a series of snapshots of the low-energy sites (LES1~5) were extracted from the simulation trajectory and compared with the initial inactive state of S1P1R. The conformation of LES5 was consistent with the three features mentioned above. Specifically, the toggle-switch on TM6 was the residue Trp269, which rotated by almost 90 degrees. The “ionic lock” between the guanidyl of residue Arg78 and the carboxyl of Glu250 was broken, and the outward displacements of TM5 and TM6 measured by the backbone carbon atom of Arg231 (TM5) and Ala245 (TM6) respectively were approximately 5.5 and 7.9 Å (Figure 4C, 4D). Therefore, this conformation was extracted as the potential active-like S1P1R model.

12 ACS Paragon Plus Environment

Page 13 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 4. (A) One-dimensional (1D) potential mean force (PMF) calculation profile for 900 ns GaMD. (B) Two-dimensional (2D) PMF calculation profile for 900 ns GaMD. (C) Alignment of conformations in the low-energy sites (the inactive S1P1R is in yellow, and the LES5 conformation is in magenta, the other conformations, LES2, LES3, and LES4, are in white). (D) Rotation of Trp269 of TM6 as the toggle-switch for activation (the inactive S1P1R is in yellow and the LES5 conformation is in magenta). In addition, the other microscopic changes were found; the Tyr221 residue on TM5 rotated inward and formed a π-alkyl interaction with the Ala139 residue on TM3. Since TM5 and TM6 shifted outward, the entrance to the intracellular domain of S1P1R was more open, leading to the subsequent entry of several water molecules. Therefore, additional hydrogen bounds were formed in water bridges between residues Tyr221 and 13 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 35

Leu135, between Asn63 and Asp91, between Asp91 and Asn303, and between Val261 and Asn307, and these bridges enhanced the interactions between the main chains of S1P1R and stabilized the intracellular domain of the receptor (Figure 5A). Moreover, the Ramachandran plot (Figure 5B) was prepared to evaluate the quality of the stabilized structure. The chosen structure obviously displayed good model quality, with 93.6% of the residues in the allowed region, 5.8% in the marginal region and only 0.6% of the residues in the disallowed region.

14 ACS Paragon Plus Environment

Page 15 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 5. (A) The water bridges inside the active-like state of S1P1R. (B) Ramachandran plot of the stabilized active-like S1P1R model (the allowed region is within the cyan curve, and the marginal region is within the magenta curve). Verification of the active-like S1P1R model by docking with representative S1P1R agonists. Two representative agonists of S1P1R, compounds 1 and 2, were examined in silico and compared with available experimental mutagenesis data. The docking results showed that the guanidyl of Arg120 and the carboxyl of Glu121 formed a salt bridge with the phosphate group and an amino group of compound 1, respectively (Figure 6A). Moreover, the aliphatic chain and phenyl ring of compound 1 were accommodated within the hydrophobic cavity of S1P1R, forming lipophilic interactions with side chains of Met124, Phe125, Trp269, Leu272, and Leu297, and these residues were key contributors to ligand binding. These docking results were consistent with the sitemutagenesis results, indicating the predicted binding mode is reliable.38-40 On the other hand, in the docking pose of oxadiazole 2 (Figure 6B), the proximal disubstituted benzene ring was adjacent and parallel to the benzene ring of the Phe125 residue, allowing them to form π-π stacking interactions. These π-π stacking interactions were essential for the potency of the agonist, which was elucidated by a site-mutagenesis study.40 Moreover, the hydrophobic substituent (trifluoromethyl or 1trifluoromethyl ethyl) and the oxadiazole ring of compound 2 strongly interacted with residues Met124, Phe210, Trp269, Leu272, and Leu297 by hydrophobic interactions, which was also indicated by the site-mutagenesis study of S1P1R.38-40 15 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 35

Figure 6. (A) The predicted binding mode of compound 1. (B) The predicted binding mode of compound 2. The salt bridge and π-π stacking between the ligand and binding site residues are depicted with red dashed lines. Exploration of the predicted binding mode using molecular dynamics simulations and MM-GBSA calculations With a rational, active-like S1P1R conformation available, conventional molecular dynamics simulations of S1P1R bound to its three clinical drug candidates (compounds 4, 5, and 6 in Figure 1) were carried out in combination with free energy calculations

16 ACS Paragon Plus Environment

Page 17 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

to fully understand the binding mode and key interactions of active-like S1P1R with its agonists. Although these compounds have completely different structures, they share the same key feature, namely, a proximal disubstituted aromatic ring with at least one electron-deficient group as a hydrophobic tail. The agonists were docked into the active-like state of S1P1R and then subjected to the 40 ns cMD simulations. The convergence and stability of the simulations were monitored by examining the RMSD values of the backbone atoms with respect to the initial structures. As shown in the plots (Figure 7), the structure of S1P1R exhibited good stability after 30 ns of simulation in all three cases (the RMSD values fluctuated from 2 Å to 3.5 Å). To ensure that the system was stable and well equilibrated, only the MD trajectories in the last 10 ns of the simulation were used for MM-GBSA analysis.

Figure 7. RMSD values of the backbone atoms of the active-like S1P1R model binding with its three agonists.

17 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 35

Herein, we elucidated the binding energy and the main driving force for the binding of the three agonists by comparing them using the MM-GBSA method. As shown in Table 1, the decomposition of the binding free energy indicated that for all three agonists the residues Met124, Phe125, Leu195, Cys206, Phe210, Trp269, Leu272, Leu276, and Leu297 were the most vital sites, which contributed to the binding mode of the active-like receptor. Among these residues, van der Waals interactions were the major contributors to the binding energy, especially interactions with Phe125. In addition to the residues set mentioned above, compound 4 bound to active-like S1P1R through salt bridges with Arg120 and Glu121 and van der Waals interactions with residues Try98, Leu128, Ser129, Leu298, Val301, and Ala300 (Figure 8A). Instead of interacting with Glu121, compound 5 interacted with the Arg120 residue through electrostatic interactions. 5 also required Leu128, Val194, Ala293, and Glu294 as additional sites of interaction (Figure 8B). Compound 6 preferentially interacted with Glu121 over Arg120 to form a salt bridge, and it also interacted with Val194, Leu125, Leu273, and Glu294 (Figure 8C). The results showed that although the three agonists possessed different structures, they formed similar interactions with residues Arg120 or Glu121 and with Met124, Phe125, Leu195, Cys206, Leu272, Leu276, Phe210, Trp269 and Leu297. The negative center and the positive center of the ligand formed a salt bridge with the guanidyl moiety of Arg120 and the carboxyl moiety of Glu121; the electron-deficient aromatic ring of the ligand formed π-π stacking interactions with the phenyl ring of Phe125, which was favorable for ligand and receptor binding in all three cases; additionally, the two proximal substituents on the tail of the ligand were surround 18 ACS Paragon Plus Environment

Page 19 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

by Leu195, Cys206, Leu272, and Leu276 and by Phe210 and Trp269, respectively. As shown in Figure 8D, these residues could be divided into four groups of binding site residues, and these residues were mainly responsible for the effective binding energies. The hydrophobic group (P2), including Met124, Phe125 and Leu297, was considered the main factor to the effective binding energy because of the π-π stacking interactions. The P1 group, containing Arg120 or Glu121, was also important for binding via the formation of salt bridges. Notably, two hydrophobic subpocket groups, P3, which includes Leu195, Cys206, Leu272, and Leu276, and P4, which includes Phe210 and Trp269, were both vital to the binding energy due to their cumulative lipophilic interactions, which have rarely been systematically researched in the binding mode of S1P1R. The four groups of binding site residues mentioned above were consistent with the previous site-mutagenesis study and could reveal the concrete details of the binding pocket in the active-like S1P1R model.38-40 Table 1. Energy contributions of S1P1R residues to the binding of ligands.a

Compound 4 Res.

V

Tyr98

Compound 5

E

P

T

V

P

T

Va

Ea

Pa

Na

Ta

-1.27

0.15

0.16

-0.16

-1.13

-0.09

-0.34

0.37

0.00

-0.06

-0.15

0.35

-0.31

-0.01

-0.12

Arg120

-0.46

-12.91

Glu121

-0.71

-14.88

8.92

-0.10

11.80

-0.14

-4.54

-0.33

-25.39

-3.94

-1.66

27.53

24.48

-0.05

-1.29

-0.56

26.40

-25.35

-0.07

0.41

-24.69

-0.21

0.98

-0.98

-41.09

38.90

-0.27

-3.45

Met124

-1.53

Phe125

-3.27

0.21

-0.82

0.29

-0.62

-0.15

-2.29

-1.55

-0.33

-3.93

-4.38

-1.07

0.71

-0.10

-2.02

-1.55

-0.48

0.57

-0.17

-1.63

-0.42

0.65

-0.33

-4.48

-4.10

-0.56

-0.04

-0.27

-4.96

Leu128

-1.44

-0.02

-0.23

-0.13

-1.82

-1.16

-0.91

0.78

-0.04

-1.32

-1.33

1.32

0.05

-0.11

-0.07

Ser129

-0.31

Val194

-0.24

-0.04

-0.06

-0.01

-0.72

0.84

-0.02

-0.42

-1.28

-0.79

0.32

-0.10

-1.85

-0.74

0.56

-0.47

-0.08

-0.72

-0.14

-1.04

2.06

-1.86

-0.13

-0.96

-1.05

-3.91

4.09

-0.10

-0.97

Leu195 Ile203

-1.28

0.35

-0.24

-0.90

0.05

-0.07

-0.11

-1.28

-2.14

-1.23

0.82

-0.23

-2.79

-2.20

0.77

-0.58

-0.18

-2.19

-0.04

-0.95

-0.25

0.10

-0.20

-0.01

-0.35

-0.81

1.15

-0.82

-0.01

-0.49

Cys206

-0.88

0.36

Phe210

-1.13

0.26

-0.01

-0.07

-0.60

-1.50

-0.75

0.54

-0.10

-1.81

-1.17

-1.75

0.61

-0.10

-2.42

-0.02

-0.16

-1.06

-2.11

-0.71

0.36

-0.21

-2.67

-0.87

0.11

-0.38

-0.07

-1.22

Trp269

-0.99

Leu272

-1.18

0.25

-0.07

-0.08

-0.90

-1.22

0.54

-0.50

-0.12

-1.31

-0.48

-0.65

0.38

-0.04

-0.79

0.03

-0.09

-0.14

-1.38

-0.80

-0.09

-0.05

-0.07

-1.01

-1.25

0.12

-0.07

-0.10

-1.31

a

a

N

Compound 6

E

a

a

a

a

a

a

N

a

a

19 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 35

Phe273

-0.67

-0.15

0.20

-0.09

-0.71

-0.30

0.08

-0.11

-0.01

-0.33

-1.40

-0.23

-0.12

-0.12

-1.88

Leu276

-0.72

0.09

-0.31

-0.08

-1.01

-0.65

-0.21

-0.01

-0.11

-0.99

-1.08

0.10

-0.32

-0.08

-1.38

Ala293

-0.17

0.05

0.01

0.00

-0.11

-0.92

-0.75

0.81

-0.09

-0.95

-0.81

0.75

-0.10

-0.09

-0.26

Glu294

-0.72

-2.99

3.19

-0.10

-0.63

-0.94

21.48

-21.44

-0.14

-1.03

-0.71

-21.36

20.96

-0.15

-1.26

Leu297

-3.55

-0.03

0.66

-0.36

-3.28

-2.54

-1.55

1.28

-0.32

-3.13

-2.31

1.70

-0.75

-0.27

-1.63

Val298

-0.90

0.07

-0.19

-0.08

-1.11

-0.10

-0.36

0.37

0.00

-0.10

-0.07

0.36

-0.39

0.00

-0.10

Ala300

-0.83

-0.10

-0.01

-0.05

-0.98

-0.11

-0.73

0.60

0.00

-0.23

-0.07

0.72

-0.72

0.00

-0.07

Val301

-0.92

0.12

-0.44

-0.07

-1.31

-0.12

-0.64

0.46

0.00

-0.30

-0.09

0.58

-0.60

0.00

-0.11

aAll

values are given in kcal/mol. The column names V: van der Waals energy contributions, E: electrostatic energy contributions, P: polar solvation energy contributions, N: nonpolar solvation energy contributions, and T: total free energy contributions.

20 ACS Paragon Plus Environment

Page 21 of 35

Journal of Chemical Information and Modeling

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

21 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 35

Figure 8. (A) The predicted binding mode of compound 4 after cMD simulation. (B) The predicted binding mode of compound 5 after cMD simulation. (C) The predicted binding mode of compound 6 after cMD simulation. (D) Per-residue contributions to the effective binding energy of the ligands, which were calculated by the MM-GBSA decomposition method. The data for compounds 4, 5, and 6 are shown in cyan, yellow and red, respectively. In addition, the P1~P4 were the residues clusters around the binding pocket. The validation of the active-like S1P1R model based on the computational and experimental approaches Three compounds containing oxadiazole linker (14a-14c) were selected from Liu et al’s manuscript,29 and modified on the polar head and the hydrophobic tail(14d, 13e) to validate the four groups of binding site residues. As depicted in Table 2, all the compounds were evaluated for activity against S1P1R in an HTRF-IP1 functional assay.41,

42

The compounds with flexible polar headgroups and electron-deficient

aromatic rings were predicted to possess the strongest interactions with S1P1R because of the formation of salt bridges with both the guanidyl moiety of Arg120 and the carboxyl moiety of Glu121 and the enhancement of the π-π stacking interactions between the phenyl ring of Phe125 and the electron-deficient aromatic ring on the tail. As expected, compound 14a, with a polar headgroup and an electron-deficient aromatic ring, not only formed salt bridge and π-π stacking interactions with the receptor but also formed hydrophobic interactions with subpockets P3 and P4 through its trifluoromethyl

22 ACS Paragon Plus Environment

Page 23 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

and isopropoxy groups. Therefore, of the five compounds, 14a exhibited the best agonist efficiency, with an EC50 value of 17.1 nM. Compound 13e, in which a propyl group was used in place of the carboxylic acid groups of compound 14a, was synthesized to verify the salt bridge interaction with P1. This modification prevented formation of a salt bridge with residue Arg120 and led to a remarkable decrease in the potency with an EC50 value of 2417 nM. In addition, compound 14b was synthesized to probe the π-π stacking interactions with Phe125 in group P2. The trifluoromethyl moiety in compound 14a was replaced by a methyl group in 14b, which preserved the hydrophobic interactions but slightly weakened the π-π stacking interactions. Therefore, the agonist activity of compound 14b was comparable to that of compound 14a, with an EC50 value of 84.7 nM. In addition, compound 14d, with a methoxy substituent instead of an isopropoxy substituent, was synthesized to explore the effect of minimizing the hydrophobic interactions with P4, and this compound showed an approximately 40-fold decrease in potency compared with compound 14a with an EC50 value of 646.5 nM. However, for compound 14c, the removal of the trifluoromethyl groups not only prevented the formation of hydrophobic interactions with group P3 but also disfavored the π-π stacking interactions with group P2, resulting in a >150-fold decrease in the potency and an EC50 value of 2762 nM. Table 2. The structures of the five oxadiazole derivatives and the results of the biological assay.a No.

Molecular structures

Glide Gscore

EC50 (nM) for S1P1R

23 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

F3 C

14a

O N

O

O N

14b

O

O

F3 C

14d

O

13e a

O

O N N

H N

COOH

H N

COOH

H N

COOH

Cl

N

F3 C

COOH

-8.425

17.1

-8.148

84.7

-7.246

2762

-8.031

646.5

-7.544

2417

Cl

N

O N

H N

Cl

N

O N

14c

Cl

N

Page 24 of 35

Cl H N

EC50 is the mean of three experimental determinations.

Considering both the docking results and the biological evaluation data, the salt bridge between the polar head group with P1, the π-π stacking interaction with the electron-deficient aromatic ring and P2, and the hydrophobic interaction between the tail of the ligand and groups P3 and P4 were identified as the four crucial factors affecting the potency of the S1P1R agonists (Figure 9A). In the further study, 1145 S1P1R agonists collected from BindingDB were docked into the inactive crystal structure of S1P1R and the active-like state of S1P1R respectively.43 As showed in Figure 9B, the correlation coefficient between the docking score and the experimental pEC50 for the agonists in the active-like state was obviously higher than in the inactive state of S1P1R. It was indicated that the active-like S1P1R model was more rational to guild the design of the S1P1R agonists.

24 ACS Paragon Plus Environment

Page 25 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 9. (A) The predicted binding mode between the five synthesized compounds and the four groups of binding site residues, P1, P2, P3, and P4 (P1 is in green, P2 is in orange, P3 is in blue, and P4 is in magenta). (B) The correlation between the experimental pEC50 values and the Glide Gscores of the 1145 agonists docked into the inactive state and active-like state S1P1R respectively. CONCLUSION In this paper, the active-like state of S1P1R was determined via GaMD simulations of the S1P1R-compound 3 complex. Then, the active-like S1P1R model was verified by docking two representative S1P1R agonists, compounds 1 and 2, into the active pocket 25 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 35

of the receptor. With a reliable active-like state of S1P1R in hand, the details of the binding pocket were revealed, and four groups of binding site residues were found using molecular dynamics simulations and MM-GBSA method for three S1P1R agonists, compounds 4, 5, and 6, in clinical trials. Group P1 included Arg120 and Glu121, which could form salt bridges with the polar headgroup of the ligand; group P2, including Met124, Phe125 and Leu297, not only formed hydrophobic interactions with the receptor but also formed π-π stacking interactions between residue Phe125 and the aromatic ring of the tail of the ligand, which was the most important interaction in the binding mode; group P3, with Leu195, Cys206, Leu272 and Leu276, was one of the hydrophobic subpockets, and group P4, with Phe210 and Trp269, was the other hydrophobic subpocket. The results of the predicted binding mode study with the four groups of binding site residues were verified by relationship between docking scores and in vitro S1P1R activities of five structural analogues synthesized in-house and 1145 known S1P1R agonists.

ASSOCIATED CONTENT Supporting Information. The Supporting Information is available free of charge on the ACS Publications website. Supporting figures for additional results of the GaMD simulations and molecular docking studies (Figures S1~3. PDF). 26 ACS Paragon Plus Environment

Page 27 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

1H

NMR, 13C NMR and HRMS spectral information for representative compounds

and biological evaluation. (PDF) AUTHOR INFORMATION Corresponding Author *E-mail: [email protected]. ORCID Xiaojian Wang: 0000-0002-1856-8820 Acknowledgements This work was financially supported by National Natural Science Foundation of China (NSFC no. 81473096 and 81602949) and CAMS Collaborative Innovation Project (no. 2017-I2M-3-011). We thank Dr. Yadong Chen for the support of the docking studies. Notes The authors declare no competing financial interest. REFERENCES (1) Roberts, E.; Guerrero, M.; Urbano, M.; Rosen, H. Sphingosine 1-phosphate Receptor Agonists: A Patent Review (2010–2012). Expert Opin. Ther. Pat. 2013, 23, 817-841. (2) Buzard, D. J.; Thatte, J.; Lerner, M.; Edwards, J.; Jones, R. M. Recent Progress in the Development of Selective S1P1 Receptor Agonists for the Treatment of Inflammatory and Autoimmune Disorders. Expert Opin. Ther. Pat. 2008, 18, 11411159. 27 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 35

(3) D'Ambrosio, D.; Freedman, M. S.; Prinz, J. Ponesimod, A Selective S1P1 Receptor Modulator: A Potential Treatment for Multiple Sclerosis and Other Immunemediated Diseases. Ther. Adv. Chronic Dis. 2016, 7, 18-33. (4) Rosen, H.; Stevens, R. C.; Hanson, M.; Roberts, E.; Oldstone, M. B. Sphingosine1-phosphate and Its Receptors: Structure, Signaling, and Influence. Annu. Rev. Biochem. 2013, 82, 637-662. (5) Pham, T. H. M.; Okada, T.; Matloubian, M.; Lo, C. G.; Cyster, J. G. S1P1 Receptor Signaling Overrides Retention Mediated by Gαi-coupled Receptors to Promote T Cell Egress. Immunity. 2008, 28, 122-133. (6) Kappos, L.; Antel, J.; Comi, G.; Montalban, X.; Oconnor, P.; Polman, C. H.; Haas, T.; Korn, A. A.; Karlsson, G.; Radue, E. W. Oral Fingolimod (FTY720) for Relapsing Multiple Sclerosis. N. Engl. J. Med. 2006, 355, 1124-1140. (7) Gonzalez-Cabrera, P. J.; Jo, E.; Sanna, M. G.; Brown, S.; Leaf, N.; Marsolais, D.; Schaeffer, M. T.; Chapman, J.; Cameron, M.; Guerrero, M. Full Pharmacological Efficacy of A Novel S1P1 Agonist That Does Not Require S1P-like Headgroup Interactions. Mol. Pharmacol. 2008, 74, 1308-1318. (8) Yamamoto, R.; Okada, Y.; Hirose, J.; Koshika, T.; Kawato, Y.; Maeda, M.; Saito, R.; Hattori, K.; Harada, H.; Nagasaka, Y., ASP4058, A Novel Agonist for Sphingosine 1-phosphate Receptors 1 and 5, Ameliorates Rodent Experimental Autoimmune Encephalomyelitis with A Favorable Safety Profile. PLoS One. 2014, 9, e110819.

28 ACS Paragon Plus Environment

Page 29 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

(9) Buzard, D. J.; Kim, S. H.; Lopez, L.; Kawasaki, A.; Zhu, X.; Moody, J.; Thoresen, L.; Calderon, I.; Ullman, B.; Han, S. Discovery of APD334: Design of A Clinical Stage Functional Antagonist of the Sphingosine-1-phosphate-1 Receptor. ACS Med. Chem. Lett. 2014, 5, 1313-1317. (10) Urbano, M.; Guerrero, M.; Rosen, H.; Roberts, E. Modulators of the Sphingosine 1-phosphate Receptor 1. Bioorg. Med. Chem. Lett. 2013, 23, 6377-6389. (11) Hanson, M. A.; Roth, C. B.; Jo, E.; Griffith, M. T.; Scott, F.; Reinhart, G.; Desale, H.; Clemons, B.; Cahalan, S. M.; Schuerer, S. C. Crystal Structure of A Lipid G Protein-coupled Receptor. Science. 2012, 335, 851-855. (12) Miao, Y.; Feher, V. A.; Mccammon, J. A. Gaussian Accelerated Molecular Dynamics: Unconstrained Enhanced Sampling and Free Energy Calculation. J. Chem. Theory Comput. 2015, 11, 3584-3595. (13) Pan, S.; Gray, N. S.; Gao, W.; Mi, Y.; Fan, Y.; Wang, X.; Tuntland, T.; Che, J.; Lefebvre, S.; Chen, Y. Discovery of BAF312 (Siponimod), A Potent and Selective S1P Receptor Modulator. ACS Med. Chem. Lett. 2013, 4, 333-337. (14) Xu, J.; Gray, F.; Henderson, A.; Hicks, K.; Yang, J.; Thompson, P.; Oliver, J. Safety, Pharmacokinetics, Pharmacodynamics, and Bioavailability of GSK2018682, A Sphingosine-1-phosphate Receptor Modulator, in Healthy Volunteers. Clin. Pharm. Drug. Dev. 2014, 3, 170-178.

29 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 35

(15) Scott, F.; Clemons, B.; Brooks, J.; Brahmachary, E.; Powell, R.; Dedman, H.; Desale, H.; Timony, G.; Martinborough, E.; Rosen, H. Ozanimod (RPC1063) is A Potent Sphingosine-1-phosphate Receptor-1(S1P1) and Receptor-5(S1P5) Agonist with Autoimmune Disease-modifying Activity. Br. J. Pharmacol. 2016, 173, 1778-1792. (16) Gohlke, H.; Case, D. A. Converging Free Energy Estimates: MM-PB(GB)SA Studies on the Protein-protein Complex Ras-Raf. J. Comput. Chem. 2004, 25, 238-50. (17) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L. T.; Lee, M. R.; Lee, T.; Duan, Y.; Wang, W. Calculating Structures and Free Energies of Complex Molecules: Combining Molecular Mechanics and Continuum Models. Acc. Chem. Res. 2000, 33, 889-897. (18) Bairoch, A.; Apweiler, R.; Wu, C. H.; Barker, W. C.; Boeckmann, B.; Ferro, S.; Gasteiger, E.; Huang, H.; Lopez, R.; Magrane, M. The Universal Protein Resource (UniProt). Nucleic Acids Res. 2005, 33, D154-D159. (19) Friesner, R. A.; Banks, J. L.; Murphy, R. B.; Halgren, T. A.; Klicic, J. J.; Mainz, D. T.; Repasky, M. P.; Knoll, E. H.; Shelley, M.; Perry, J. K. Glide: A New Approach for Rapid, Accurate Docking and Scoring. 1. Method and Assessment of Docking Accuracy. J. Med. Chem. 2004, 47, 1739-1749. (20) Jo, S.; Klauda, J. B.; Im, W. CHARMM-GUI Membrane Builder for Mixed Bilayers and Its Application to Yeast Membranes. Biophys J. 2009, 97, 50-58.

30 ACS Paragon Plus Environment

Page 31 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

(21) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696-3713. (22) Dickson, C.; Madej, B. D.; Skjevik, Å. A.; Betz, R. M.; Teigen, K.; Gould, I. R.; Walker, R. C. Lipid14: The Amber Lipid Force Field. J. Chem. Theory Comput. 2014, 10, 865-879. (23) Miao, Y.; Feher, V. A.; Mccammon, J. A. Gaussian Accelerated Molecular Dynamics: Unconstrained Enhanced Sampling and Free Energy Calculation. J. Chem. Theory Comput. 2015, 11, 3584-3595. (24) Pang, Y. T.; Miao, Y.; Wang, Y.; Mccammon, J. A. Gaussian Accelerated Molecular Dynamics in NAMD. J. Chem. Theory Comput. 2017, 13, 9-19. (25) Miao, Y.; Huang, Y. M.; Walker, R. C.; Mccammon, J. A.; Chang, C. A. Ligand Binding Pathways and Conformational Transitions of the HIV Protease. Biochemistry. 2018, 57, 1533-1541. (26) Roe, D. R.; Rd, C. T. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9, 3084-3095. (27) Miao, Y.; Sinko, W.; Pierce, L.; Bucher, D.; Walker, R. C.; Mccammon, J. A. Improved Reweighting of Accelerated Molecular Dynamics Simulations for Free Energy Calculation. J. Chem. Theory Comput. 2014, 10, 2677-2689. 31 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(28) Fulle, S.; Withers-Martinez, C.; Blackman, M. J.; Morris, G. M.; Finn, P. W. Molecular Determinants of Binding to the Plasmodium Subtilisin-like Protease 1. J. Chem. Inf. Model. 2013, 53, 573-583. (29) Liu, T.; Jin, J.; Chen, Y.; Xi, Q.; Hu, J.; Jia, W.; Chen, X.; Li, Y.; Wang, X.; Yin, D. Identification and Structure–Activity Relationship (SAR) of Potent and Selective Oxadiazole-based Agonists of Sphingosine-1-phosphate Receptor (S1P1). Bioorg. Chem. 2019, 82, 41-57. The research work about the induction of the active-like state of S1P1R mentioned in Liu’s article was not described yet. which is done in the current manuscript. (30) Liao, J. M.; Wang, Y. T. In Silico Studies of Conformational Dynamics of Mu Opioid Receptor Performed Using Gaussian Accelerated Molecular Dynamics. J. Biomol. Struct. Dyn. 2018, 1-23. (31) Schwartz, T. W.; Frimurer, T. M.; Holst, B.; Rosenkilde, M. M.; Elling, C. E. Molecular Mechanism of 7TM Receptor Activation--A Alobal Toggle Switch Model. Annu. Rev. Pharmacol. Toxicol. 2006, 46, 481-519. (32) Holst, B.; Nygaard, R.; Valentinhansen, L.; Bach, A.; Engelstoft, M. S.; Petersen, P. S.; Frimurer, T. M.; Schwartz, T. W. A Conserved Aromatic Lock for The Tryptophan Rotameric Switch in TM-VI of Seven-transmembrane Receptors. J. Biol. Chem. 2010, 285, 3973-3985. (33) Ballesteros, J. A.; Jensen, A. D.; Liapakis, G.; Rasmussen, S. G.; Shi, L.; Gether, U.; Javitch, J. A. Activation of the Beta 2-adrenergic Receptor Involves Disruption of 32 ACS Paragon Plus Environment

Page 32 of 35

Page 33 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

An Ionic Lock Between the Cytoplasmic Ends of Transmembrane Segments 3 and 6. J. Biol. Chem. 2001, 276, 29171-29177. (34) Rosenbaum, D. M.; Rasmussen, S. G.; Kobilka, B. K. The Structure and Function of G-protein-coupled Receptors. Nature. 2009, 459, 356-63. (35) Deupi, X.; Standfuss, J. Structural Insights into Agonist-induced Activation of G-protein-coupled Receptors. Curr. Opin. Struct. Biol. 2011, 21, 541-551. (36) Scheerer, P.; Park, J. H.; Hildebrand, P. W.; Yong, J. K.; Krau, N.; Choe, H. W.; Hofmann, K. P.; Ernst, O. P. Crystal Structure of Opsin in Its G-protein-interacting Conformation. Nature. 2008, 455, 497-502. (37) Park, J. H.; Scheerer, P.; Hofmann, K. P.; Choe, H.; Ernst, O. P. Crystal Structure of The Ligand-free G-protein-coupled Receptor Opsin. Nature. 2008, 454, 183-187. (38) Fujiwara, Y.; Osborne, D. A.; Walker, M. D.; Wang, D. A.; Bautista, D. A.; Liliom, K.; Van Brocklyn, J. R.; Parrill, A. L.; Tigyi, G. Identification of The Hydrophobic Ligand Binding Pocket of the S1P1 Receptor. J. Biol. Chem. 2007, 282, 2374-2385. (39) Naor, M. M.; Walker, M. D.; Van Brocklyn, J. R.; Tigyi, G.; Parrill, A. L. Sphingosine 1-phosphate pKa and Binding Constants: Intramolecular and Intermolecular Influences. J. Mol. Graphics Model. 2007, 26, 519-528.

33 ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 35

(40) van Loenen, P. B.; de Graaf, C.; Verzijl, D.; Leurs, R.; Rognan, D.; Peters, S. L. M.; Alewijnse, A. E. Agonist-dependent Effects of Mutations in The Sphingosine-1phosphate Type 1 Receptor. Eur. J. Pharmacol. 2011, 667, 105-112. (41) Norskovlauritsen, L.; Thomsen, A. R.; Braunerosborne, H. G Protein-Coupled Receptor Signaling Analysis Using Homogenous Time-Resolved Förster Resonance Energy Transfer (HTRF®) Technology. Int. J. Mol. Sci. 2014, 15(2), 2554-2572. (42) Degorce, F.; Card, A.; Soh, S.; Trinquet, E.; Knapik, G. P.; Xie, B. HTRF: A Technology Tailored for Drug Discovery-A Review of Theoretical Aspects and Recent Applications. Current Chemical Genomics. 2009, 3(1), 22-32. (43) Liu, T.; Lin, Y.; Wen, X. BindingDB: A Web-accessible Database of Experimentally Determined Protein-ligand Binding Affinities. Nucleic Acids Res. 2007, 198-201.

34 ACS Paragon Plus Environment

Page 35 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

For Table of Contents Use Only

A Computational Approach to the Study of the Binding Mode of S1P1R Agonists Based on the Active-like Receptor Model Yonghui Chena, Tianqi Liub, Qiumu Xib, Wenqiang Jiaa, Dali Yina,b, Xiaojian Wang*a,b

35 ACS Paragon Plus Environment