Rational Design of Highly Potent and Slow-Binding Cytochrome bc1 Inhibitor as Fungicide by Computational Substitution Optimization

Hit to lead (H2L) optimization is a key step for drug and agrochemical discovery. A critical challenge for H2L optimization is the low efficiency due to the lack of predictive method with high accuracy. We described a new computational method called Computational Substitution Optimization (CSO) that has allowed us to rapidly identify compounds with cytochrome bc1 complex inhibitory activity in the nanomolar and subnanomolar range. The comprehensively optimized candidate has proved to be a slow binding inhibitor of bc1 complex, ~73-fold more potent (Ki = 4.1 nM) than the best commercial fungicide azoxystrobin (AZ; Ki = 297.6 nM) and shows excellent in vivo fungicidal activity against downy mildew and powdery mildew disease. The excellent correlation between experimental and calculated binding free-energy shifts together with further crystallographic analysis confirmed the prediction accuracy of CSO method. To the best of our knowledge, CSO is a new computational approach to substitution-scanning mutagenesis of ligand and could be used as a general strategy of H2L optimisation in drug and agrochemical design.

Agrochemicals play a significant role in modern agriculture by protecting crops against yield loss. Due to the increasing demand for food supply and the explosive development of resistance, the discovery and development of new agrochemicals are continuously demanded 1 . However, agrochemical discovery remains a vital challenge 2 : the research and development (R&D) costs have been rising from U.S. $152 million in 1995 to $256 million in 2005 and the number of compounds needs to be synthesized have been rising from 52500 in 1995 to 140000 in 2005 to bring a new agrochemical to market. A new agrochemical is invented initially by the discovery of a bioactive lead compound. Thus, lead generation is a crucial step for the agrochemical discovery 3 , which requires at least two processes: hit identification and Hit-to-Lead (H2L) optimization 4 . Over the past decade, introduction of new technologies, such as high-throughput screening (HTS) and fragment-based design, have made hit identification become more efficient 5,6 . But, traditional H2L optimization is still carried out through a large number of analogues' syntheses, which lead to a low efficient and resource-intensive research 7 . Although some techniques, such Scientific RepoRts | 5:13471 | DOi: 10.1038/srep13471 as quantitative structure-activity relationship (QSAR) 8 and molecular dynamics (MD) simulations 9,10 , have been developed to help understand the structure-activity relationship for H2L optimization, they still have some drawbacks, such as dependent of sufficient biological data, highly time-consuming, and low prediction accuracy 11 . Therefore, there is urgent demand in developing an effective tool enabling to eliminate undesirable lead compound classes/types early during H2L optimization process, prior to expensive and extensive experimental work 12 . Having this challenge in mind, we developed a new computational protocol called Computational Substitution Optimization (CSO) to rapidly assess how replacing hydrogen atoms with new substituent into a "hit" scaffold could improve binding free energies to their targets. To showcase this method, we applied CSO in the H2L optimization of cytochrome bc 1 complex inhibitors to discover antifungal lead.
The cytochrome bc 1 complex (EC 1.10.2.2, bc 1 ) has been chosen as a model for the present study, because of its crucial role in the life cycle. The bc 1 complex is an essential respiratory enzyme complex present in the inner mitochondrial membrane of eukaryotic organisms. It has been identified as a promising target for new drugs and agrochemicals 13 . For example, malaria, as a devastating tropical disease, remains a major threat to over 40% of the world's population and approximately 660,000 deaths in 2010 14 . Atovaquone is a major drug for prophylaxis and treatment of uncomplicated malaria patient 15 and as alternative therapy in case of resistances against chloroquine or artemisin-based therapies 16 . But with a limited number of antimalarial drugs available, the emergence of multidrug-resistant Plasmodium falciparum malaria threatens the health of individual patients as well as eradication strategies for the disease 17 . In addition, some bc 1 complex inhibitors have been introduced into the agrochemical market to provide effective control of important fungal diseases threatening global food security 18 . Azoxystrobin (AZ), as a representative inhibitor targeting bc 1 , is the biggest using fungicide in the whole word with a global annual sales of over USD1.4 billion 19 . However, the increasing number of resistance is also a big problem that AZ has to face 19 . Therefore, there is a wide research interest on the discovery of cytochrome bc 1 complex inhibitor today, which could be used as a specific probe for novel functional study and as a new lead compound for future drug and agrochemical discovery 20 .
In this study, the application of CSO has led to a series of new bc 1 inhibitors with nanomolar to subnanomolar potency. Excellent agreements between the predicted and experimental binding free energies were obtained, indicating the predictive capacity of CSO. Further, assays showed that compound 18 was the lead candidate from the new series of bc 1 complex inhibitors because of its high potency (K i = 4.10 nM which is up to 73-fold higher than that of AZ), outstanding physiochemical properties, and excellent in vivo bioactivity. In order to investigate the underlying mechanisms responsible for the excellent in vivo bioactivity, the kinetic effects on cytochrome bc 1 complex were studied and exhibited slow-binding characteristics, which is different from the classical fast-binding of AZ indicating a turning point for bc 1 complex inhibitor discovery. Most importantly, the determined crystal structure of lead compound bound to bc 1 complex at 3.23 Å provided a solid molecular basis for the prediction accuracy of CSO method. This work demonstrated that CSO method can significantly improve the efficiency of H2L optimization process and could be used as a general strategy for drug and agrochemical discovery.

Results
H2L Optimization through CSO. The workflow of the computational protocol used in this study is shown in Fig. 1. This new protocol is a combination of a MD simulation and a free-energy perturbation (FEP)-based scanning calculation. Due to E-Methyl 2-(2-((benzo-thiazol-2-ylthio)methyl)phenyl)-3-methoxyacrylate (compound 1) was identified previously as a hit compound targeting cytochrome bc 1 complex (K i = 31.10 nM, against porcine bc 1 ) 21 , the MD simulation was first performed on compound 1-bound bc 1 complex system. To evaluate the convergence of the MD simulation, plots of root-meansquare deviation (RMSD) and the distances of hydrogen bonds were examined (Fig. S1).
Due to direct interactions of the substituent with the phenyl ring of phenylalanine, electron-donating and electron-withdrawing substituent can effectively increase the π -π stacking interactions 22 (Table 1). To avoid excessive perturbation in FEP-based calculation, only single substitutions were investigated in the present study. In total, 80 computational substitutions were made and the computational changes in binding free energies as a result of each substitution were calculated for each compound. Detailed results of the binding free-energy changes are provided in Table S1.
A 10-fold improvement of potency should be the minimum standard for an acceptable H2L optimization process. Thus, by setting the criterion of the calculated binding free-energy shift (∆∆G cal ) lower than − 1.37 kcal/mol (~10-fold improvement in activity), 29 compounds were considered as candidates for further chemical synthesis. Considering the feasibility of organic synthesis, we only synthesized 10 out of them which are easy to synthesize. Their experimental binding free-energy shifts (∆∆G exp ) ranged from − 2.69 to − 1.40 kcal/mol. For an ideal method, it can predict not only compounds with higher potency, but also compounds with lower potency. Hence, a further 11 compounds those ∆∆G cal larger than − 1.37 kcal/mol and easier to synthesise were also selected randomly for chemical synthesis. The synthetic routes for the total 21 compounds are summarised in Scheme S1. Their chemical structures were characterised by 1 H and 13 C NMR spectroscopy, HRMS, elemental analysis, and X-ray diffraction analysis (Fig. S2).
In Vitro Activities and Lead Identification. Inhibitory kinetics is of great importance to understand the inhibitory activity of compounds. The inhibition constants (K i = 0.39-159.96 nM) of these compounds with the bc 1 complex were determined according to previously established method 23 and listed in Table 1. These results indicate that the newly synthesized compounds are specific and effective inhibitors of bc 1 complex. Notably, the experimental and calculated binding free-energy shifts (Δ Δ G exp and Δ Δ G cal , respectively) of these compounds showed good linear correlation (r 2 = 0.8), confirming the predictive accuracy of the CSO method.
Based on the above in vitro study, compounds with nanomolar and subnanomolar K i values were selected from the total 21 newly synthezised compounds for further studies. Pesticide likeness 24 , defined as a complex balance of various physicochemical properties and structure features, can evaluate whether a molecule is similar to the known pesticides. After the filtration of pesticide likeness, seven compounds were kept to adhere to our previously defined rules for pesticide-likeness including molecular weight (MW) ≤ 435 Da, ClogP ≤ 6, number of H-bond acceptors (HBA) ≤ 6, number of H-bond donors (HBD) ≤ 2, number of rotatable bonds (ROB) ≤ 9, and number of aromatic bonds (ARB) ≤ 17 (Table 1) 24 . Finally, based on the results of in vivo test, three out of seven compounds were found to display good antifungal activity (> 70%) against downy mildew and powdery mildew at the concentration of 200 mg/L (Table S2). Among these three compounds, compound 18 completely inhibited the growth of both species of fungus, outperforming the control fungicide, AZ (which inhibited the growth of powdery mildew, but only reduced the growth of downy mildew by 94%). Therefore, compound 18 was recognized as lead compound for subsequent studies.
In Vivo Antifungal Activity. The above results indicated that compound 18 shows excellent potency towards bc 1 complex, pesticide-like properties, and antifungal activity. For a further validation, the controlling efficacy of compound 18 against downy mildew and powdery mildew disease was tested in greenhouse and field models. Because the pathogens exhibit a specific host for cucumber, the inoculations of Pseudoperonospora cubensis (P. cubensis) and Sphaerotheca fuliginea (S. fuliginea) were carried out by spraying a conidial suspension on the seed leaves of cucumber. As expected, the greenhouse test indicated that compound 18 showed excellent fungicidal activity against both downy mildew and powdery mildew diseases. The protective and curative EC 90 values of compound 18 against P. cubensis were 4.29 mg/L and 120.29 mg/L, respectively, and against S. fuliginea are 1.89 mg/L and 32.89 mg/L, (1) The MD simulation of protein with original "hit" ligand was performed. (2) In order to obtain a representative ensemble of the binding structures, snapshots were collected from the MD trajectory at regular intervals and minimized. (3) The "Hit" ligand in each snapshot collected from the MD trajectory was mutated to new ligand with substitutions using the library of substitution and minimized. Each substitution group was set to link with the specific site of "Hit" compounds. (4) Δ Δ G was calculated and the final binding free energy change is the average of the calculated values associated with each snapshot. New leads were finally determined according to the order of the binding free-energy shifts.
respectively. As a control, the protective and curative EC 90 values of AZ against P. cubensis are 5.36 mg/L and 50.66 mg/L, respectively, and against S. fuliginea are 254.30 mg/L and > 500 mg/L, respectively.
To further study the potential of compound 18 against powdery mildew and downy mildew, field experiments were conducted during the growing season of summer squash and cucumber. The summer squash plants at a fruiting stage naturally infected by S. fuliginea were used. A foliar application of compound 18 at 75 g.ai/ha recorded 74.60% disease reduction over control. Almost no lesions were observed on the leaves of summer squash. If the treatment concentration of compound 18 was reduced by half to 37.5 g.ai/ha, lesions on the leaves of treated plants expanded very slowly, or ceased to expand, which demonstrates that compound 18 has a curative activity at 37.5 g.ai/ha against the expansion of lesions by S. fuliginea with a curative effect of approximately 70.92% ( Fig. 2A,B). Lesions on the leaves of control plants, however, expanded quickly, and intense fungal hyphae could be observed. The curative effect is approximately 64.85% for the AZ-treated (93.75 g.ai/ha) summer squash. In addition, field trials for the treatment of downy mildew of cucumber were also performed. Cucumber plants naturally infected by P. cubensis were used for the study. Fungi on the leaves of control plants expanded quickly and lesions could be observed. However, we observed that when compound 18 (60 g.ai/ha) was sprayed on the leaves of the cucumber at the primary infection location, only small or no lesions were observed, and the statistical curative effect was approximately 82.59% (Fig. 2B). As a comparison, when leaves of cucumber were treated with AZ (75 g.ai/ha), a lower curative effect (59.59%) was observed. The field experiments show that the curative effect of compound 18 is superior to AZ as a treatment of powdery mildew and downy mildew.
Action Mechanism. Based on the above in vitro and in vivo studies, it is reasonable to speculate that compound 18 might have different action mechanism from AZ. Hence, the kinetic effects of compound 18 on cytochrome bc 1 complex were investigated and compared with compound 1 (hit compound) and AZ (Fig. 3A).
Our previous data showed that AZ inhibits the activity of SCR with K i = 297.6 ± 7.9 nM, and that AZ is non-competitive with respect to cytochrome c 23 . To address the competition of AZ with respect to ubiquinol, we conducted inhibition assays for bc 1    in the absence and presence of the inhibitor AZ. Figure 3B shows a set of double-reciprocal plots of 1/v 0 versus 1/[DBH 2 ], which yields a series of straight lines with a common intercept on the 1/v axis with different slopes. Because the intercepts of each line on the 1/[DBH 2 ] and 1/v axes are equal to − 1/K m and 1/V max , respectively, at the indicated AZ concentration, we can conclude that K m for the DBH 2 substrate increases with increasing AZ concentration, whereas V max is not affected. These results indicate that AZ is a competitive inhibitor of bc 1 complex with respect to substrate DBH 2 , and its kinetic parameters are determined as K m(DBH2) = 80.0 ± 7.0 μ M and K i = 394.7 ± 19.6 nM.
To unravel the inhibitory mechanism of compound 1, we first examined the dependence of product formation on the concentration of the cytochrome c as substrate by monitoring the initial rates of SCR (Assay 1). When the succinate concentration was fixed and that of cytochrome c varied in the presence of various fixed concentrations of compound 1, the double-reciprocal plots yielded a non-competitive pattern (Fig. 3C). Analysis of the data yielded K m and K i values of 4.3 ± 0.2 μ M and 31.1 ± 0.9 nM, respectively. Next, we analysed the effect of ubiquinol as substrate in the DBH 2 -cytochrome c system (Assay 3). When DBH 2 was used as the substrate, a competitive pattern was observed for compound 1, in which a series of straight lines converged to a point on the vertical axis (Fig. 3D). The K i value was determined as 96.5 ± 5.6 nM, which is approximately 3-fold higher than that determined in the succinate-cytochrome c system. This difference likely results from the presence of lauryl maltoside in the assay system, as suggested in previous studies 23 . Taken together, compound 1 is a non-competitive inhibitor with respect to substrate cytochrome c, but a competitive inhibitor with respect to DBH 2 , similar to fungicide AZ. Different with AZ and compound 1, compound 18 unexpectedly behaves as a slow-binding methoxyacrylate (MOA)-type inhibitor of bc 1 complex (Figs 4A and S3). In the presence of compound 18, the product formation curve for the SCR-catalysed reaction begins linearly but falls off and approaches a steady state with increasing time. However, when the enzyme was pre-incubated with compound 18 to reach the binding equilibrium, the progress curve simply displays a linear product-versus-time relationship (line 3), with a rate approximately the same as the steady-state rate of line 2 (Fig. S3). These data suggest that compound 18 is a slow-binding inhibitor of bc 1 complex.
To reveal the kinetic mechanism underlying this inhibition, we first examined the effect of the concentration of the cytochrome c as substrate and the inhibitor compound 18 on product formation by following the SCR-catalysed reaction (Assay 1). In the presence of different concentrations of cytochrome c and a fixed concentration of compound 18, the time-dependence of the reactions follows the substrate reaction kinetic theory, and the kinetic parameters (the initial [v 0 ] and steady-state [v s ] rate of the reaction and the observed first-order rate constant [k obs ]) can be determined by nonlinear regression analysis ( Fig. 4A). As shown in the inset, the plot of k obs against the concentration of compound 18 is a horizontal line, which clearly indicated that compound 18 is non-competitive with respect to cytochrome c 23 . Consequently, the values of the true association and dissociation rate constants k +0 and k −0 were equal to those of the apparent constants A and B, respectively, and can be determined from one set of enzymatic assays by varying the concentration of compound 18 at a fixed cytochrome c concentration (Fig. 4B). The kinetic parameters of compound 18 were determined, in which the inset shows that k obs is proportional to the inhibitor concentration with the slope and the intercept giving k +0 (= 0.00066 ± 0.00001 s −1 nM −1 ) and k −0 (= 0.00271 ± 0.00031 s −1 ), respectively 23 . The inhibition constant for compound 18 can be calculated as K i = k −0 /k +0 = 4.1 ± 0.5 nM.
Next, we assessed the relationship between compound 18 and the other substrate, ubiquinol, by using the DBH 2 -cytochrome c system (Assay 3). At a fixed concentration of compound 18 and different concentrations of DBH 2 , the progress curves exhibited similar curvilinear patterns as observed in the assays at varying cytochrome c concentration, as a result of the slow onset of inhibition (Figs  S4 and 4A). However, the k obs values decreased with increasing concentrations of DBH 2 (inset of Fig.  S4) in contrast to the observation in the inset of Fig. 4A. To achieve a comprehensive understanding, we followed the rate of product formation at various concentrations of DBH 2 and compound 18. Two representative sets of curves were shown in Fig. 4C,D, and the k obs values at the indicated DBH 2 concentrations were determined by using the procedures described above. As shown in the insets, both sets of k obs values were proportional to the concentration of compound 18, and the slopes and intercepts of these straight lines denote the apparent association and dissociation constants (A and B, respectively) at various DBH 2 concentrations 23 . The characteristic signature of decreasing A values and invariable B values with increasing concentrations of DBH 2 suggests that compound 18 is a competitive inhibitor with respect to DBH 2 (Fig. 4E,F). With the determined K m(DBH2) = 80.0 μ M, the microscopic inhibition rate constants k +0 = 0.00038 ± 0.00001 s −1 nM −1 and k −0 = 0.00509 ± 0.00009 s −1 , and inhibition constant K i = k −0 /k +0 = 13.4 ± 0.6 nM can be derived. Together, these kinetic analyses suggest that compound 18 is a slow-binding inhibitor of bc 1 complex, and is non-competitive with substrate cytochrome c but competitive with respect to DBH 2 .

Structural Basis of Cytochrome bc 1 Complex. Elucidation of the structural basis of compound
18 inhibition is helpful to further understand its action mechanism and to validate the reliability of the predicted binding modes of inhibitors. Hence, we determined the crystal structures of compound 18 bounds to chicken bc 1 complex at a resolution of 3.23 Å. The binding mode of compound 18 is clearly defined by the electron density map in Fig. S5. According to the crystal structure, compound 18 fits well in the pocket formed by side chains Phe275, Phe129, Ile147, Pro271, Glu272, Leu295, Met125, and Val299. The structural similarity between the predicted and X-ray crystal binding models was found to be RMSD = 0.635 Å (Fig. S6), further validated the prediction capability of CSO method.
A structural comparison of the binding modes of compound 18 and AZ has revealed interesting similarities and differences. The overall arrangement of the compound 18-bc 1 complex resembles the structure of the AZ-bc 1 complex with no large conformational changes of the pocket, especially the side chain of the surrounding residues (Fig. 5A,B). Compound 18 mimics the binding mode of AZ, with the methoxyacrylate and the bridging phenyl ring overlapping extensively with that of AZ, both of which bind in a sub-pocket formed by Tyr132, Glu272, Pro271, Tyr279, and Gly143 (Fig. 5C,D). However, notable differences appear in the binding mode of the side chain. The π -π interactions between Phe275 and the benzothiazole side-chain group of compound 18 are improved greatly, relative to that of AZ, because of its large co-planar structure (Fig. 5E,F), which can induce ~73-fold improvement of the in vitro activity based on the K i value of AZ.

Discussion
An effective tool to eliminate undesirable compounds prior to expensive and extensive experimental work is crucial for improving the efficiency of H2L optimization. Upon this challenge, CSO method was developed and proved to be an effective method to improve the efficiency of H2L process. Compared with The π -π interactions between Phe275 and the benzothiazole side-chain group of compound 18. (F) Distorted π -π interactions between Phe275 and the side-chain group of AZ. Stereo views of the X-ray structure of compound 18 bound to chicken bc 1 complex are available in Fig. S5. the traditional methods, the CSO method has several advantages. First, CSO method can significantly reduce the numbers of compounds, need to be synthesized. For the present study, we need to synthesize 80 compounds theoretically, but only 21 compounds were synthesized (actually, only ten compounds needed to be synthesized, and the other 11 compounds were just for testing) and 12 compounds with nanomolar to subnanomolar potency were obtained. That means, taking the present study as an example, CSO method improved the efficiency of H2L by about 75%. Second, the CSO method is independent of biological data and molecular descriptors, and has the ability to predict a more expansive chemical space. The high correlation of Δ Δ G (r 2 = 0.8 for the linear correlation) and the high correspondence of the predicted and crystallographic poses (RMSD = 0.635 Å) support the reliability and accuracy of the CSO method. On the contrary, the traditional QSAR models are highly dependent on biological data and molecular descriptors, and always used to explain the existing results and have very poor ability to predict the results outside the training sets. Third, because the CSO method is based on the assumption that the binding modes of substituted hits are similar to that of the original hit, it is possible to perform adequate conformational sampling based on the energy-minimized proteins in complex with substituted hits. Therefore, the CSO method is significantly time-saving compared with MD models.
It should be noted that the CSO method also have limitations. Clearly, this method is based on the assumption that the binding modes of substituted hits are similar to that of the original hit. This assumption is reasonable for small substituents. For a more bulky substituent, it may cause a considerable change in the binding mode and hence leads to significantly overestimate of the binding free-energy shift. Therefore, further improvement should be done to make CSO method suitable for broader application (more bulky groups).
The application of CSO method resulted in a candidate (compound 18). Compared with the commercial fungicide (AZ), the candidate has several advantages as follows: First, compound 18 has a high inhibiting potency for bc 1 complex (K i = 4.10 nM), which is up to ~73-fold more potent than AZ (K i = 297.6 nM). In the view of ligand efficiency (LE) 25 , compound 18 (LE = 0.45 kcal mol −1 per non-H atom) is prior to AZ (LE = 0.3 kcal mol −1 per non-H atom). Second, compound 18 exhibited stronger antifungal activity against downy mildew and powdery mildew than AZ. In the greenhouse test and field trial, the protective and curative effect of compound 18 against downy mildew and powdery mildew are much higher than that of AZ. Third, the inhibitory kinetics studies showed that compound 18 exhibit slow-binding characteristics, which is great different from the classical fast-binding of AZ. This significant characteristic may help to overcome the increasing resistance problem that AZ has to face. Finally, the chemical structure of compound 18 is much simple and easily to be synthesized. The successful design of a highly efficient fungicide candidate also demonstrates a promising future of CSO method. In fact, the commercial development of compound 18 in China is in progress.

Materials and Methods
Computational protocol. Based on the modification and combination of AutoGrow and the Amber 9.0 program 26,27 , the CSO protocol was designed to perform automatically computational substitution, energy minimization, and binding affinity evaluation. Depicted in Fig. 1 is the workflow of the CSO calculations.
Kinetic assays. The porcine succinate-cytochrome c reductase (SCR, the mixture of complex II and bc 1 ) was prepared essentially according to the previously reported method 28 . The enzymatic activities of SCR, complex II, and bc 1 complex were analyzed in respective reaction mixtures as reported previously 29 . Kinetic analyses for the inhibition mechanism were performed as previously described 23 .
Greenhouse Fungicidal Activity. The fungicidal activity of compound 18 against P. cubensis and S. fuliginea was evaluated according to a previously published method 30 and a potted-plant test method was adopted.
Field Trials. In the field trials, summer squash plants which is naturally infected by S. fuliginea and cucumber plants which is naturally infected by P. cubensis were used.
Crystallization and structure determination. Orthorhombic crystal of chicken bc 1 in the space group P2 1 2 1 2 1 , containing a dimer in the asymmetric unit was prepared under optimized initial crystallization conditions. The structure was validated using the online tools at the molprobity site 31 and deposited at the protein data bank with ID 4U3F.