Abstract
Keywords
Central nervous system Citalopram Donepezil Binding Target Response Surface
Introduction
Mental disorders include a wide range of common neurological and psychiatric illnesses. The nature of CNS disorders changes across the human lifespan (1, 2) and affects a huge amount of people worldwide (3). Mental and neurological disorders pose the largest health, economic and social capital burden worldwide of any disease group. Indeed, the proportionate share of the total global burden of neurologic disorders is projected to rise, highlighting an urgent need for more selective and potent drugs to treat CNS disorders (1). The process of drug development is challenging, timeconsuming, expensive, and requires tedious steps. To overcome these problems to some extent, insilico drug design approaches seem to be cost and timeefficient procedures. In this regard, valuable information about certain CNS targets and their interacted ligands/drugs is achievable via relevant online databases. On the basis of such data sources, molecular modeling studies aiming at structural elucidation of ligandreceptor interactions could be well established and developed. The final goal of such structurebased insilico studies would be to accurately and precisely predict the interactions of candidate small molecules within desired target binding sites. According to this, molecular docking is the most popular virtual structurebased method since it predicts spatial pose(s) of ligands inside the binding site of docked receptors while at the same time estimates the affinity toward the macromolecular target(s) (4). Obtained information from the docking technique is useful for attaining drugmacromolecule complexes with optimized conformations and less binding free energy (5). Search algorithm and scoring function are principal components of different docking methods. Commonly utilized search algorithms are genetic algorithm, Monte Carlo, fragmentbased and molecular dynamics within a popular docking package such as DOCK (6), AutoDock (7), Gold (8), FlexX (9), and Glide (10). There are different approaches to molecular docking procedures, which are mainly categorized as rigid ligand/rigid target, flexible ligand/rigid target, and flexible ligand/rigid target. However, a few docking packages provide flexible ligand/flexible target via allowing some kinds of protein flexibilities such as sidechain movements.
Within recent advancements and regarding the important role of computational modeling in drug design, more correlated insilico invitro/invivo experimental data has a determinant significance. More valid insilico data with higher confidence intervals might be envisaged through careful inspection/optimization of the effective factors and their interactions on the final response. Traditional optimization methods consider the variation of one factor while holding others constant. It has been revealed that such customary techniques require more trial runs and are exclusively focused on the effect of just varied factors. This is a major technical bottleneck since the interferences among factors are not taken into account. Moreover, when multiple methodological factors are involved in a typical procedure, this technique becomes unproductive and timeconsuming. Response surface methods (RSM) have been established to study factors bearing more than three levels in which different models can be developed (11, 12). Briefly speaking, RSMs offer two distinctive advantages; simultaneous exploration of factor effects enabling to record interactive effects and requirement for lower trial runs to optimize the process. The latter issue may be very beneficial in saving time and money. Factor levels might be selected upon previous knowledge of the logic numerical or categorical range. Moreover, a quantifiable response is the most important prerequisite to run such statistical designs (design of experiments; DOE).
To our best knowledge, no reports on the application of RSM for evaluating effective parameters on molecular docking accuracy have been reported yet. Our previous work focused on the application of a multifactor RSM analysis to model a docking of fluconazole against various CYP51 conformations with the aim of identifying and ranking significant and interactive effects of computational factors on the docking output of a potent antifungal drug (13). In continuation to our interest in the relevant field (13), we aimed at unveiling the full potential of RSM approaches in optimizing molecular docking simulations with a particular focus of the current study on improving AutoDock4.2 driven binding affinities of citalopramserotonin transporter (SERT) and donepezilacetyl cholinesterase (AChE) complexes as prototype systems.
To address the rationale behind selecting citalopram and donepezil as candidate CNS drugs in this study, a few words are said here regarding the pathophysiology of relevant disorders and the clinical importance of drugs.
Alzheimer’s disease (AD) is the most common neurodegenerative disorder and the sixth most common cause of death in featuring gradually progressive cognitive and functional deficits as well as behavioral changes and is associated with accumulation of amyloid and tau depositions in the brain (14). The cholinergic nervous system and acetylcholinesterase activity are closely related to the pathogenesis of AD. The most used therapy of AD is based on enhancing cholinergic function using inhibitors of acetyl AChE like rivastigmine, donepezil, or galantamine (15). Donepezil is an important oral medication used to improve cognition and behavior in people involved with AD. Depression is a familial mood disorder that causes a persistent feeling of sadness and loss of interest. Also called Major Depressive Disorder (MDD), it affects how you feel, think and behave and can lead to a variety of emotional and physical problems. As it has been predicted that MDD would be the second leading cause of death and disability by the year 2020, it became an ideal target for pharmacogenetic approaches. Among all choices of MDD treatment, the selective serotonin reuptake inhibitor (SSRI) antidepressants such as citalopram are mentioned as the firstline treatment of depression. Genetic variation of SERT is involved in the clinical remission of major depressive episodes after citalopram treatment (16). From the pharmacological aspect of view, candidate targets were selected with regard to the most studied molecular pathways within major disorders worldwide, namely serotonin reuptake inhibition by citalopram (17) and AChE inhibition by donepezil (18) (Table1).
Reference invitro binding data was retrieved from Binding MOAD, PDB Bind, and Binding DB data banks. For this purpose, candidate drugs were docked inside the binding sites of SERT and AChE according to the BoxBehnken designed matrix. Within the assembled drug matrices, response changes were monitored with simultaneous variations of factor levels, and results were subjected to statistical analysis to produce quadratic models (Figure 1). The final step included numerical optimization with ideal docking accuracies with the aim of achieving enhanced methodological conditions with regard to financial and time restrictions. As it is obvious from the above explanations, the major aim of the current work was to study the effectiveness and accuracy of docking results for a dependent drugtarget system and hence not the comparison between drugs. On the basis of this foundation, drugs were not necessarily needed to be chosen from one category since we aimed at numerical optimization of different drugtarget interaction systems.
Experimental
Drug/target
Citalopram and donepezil were selected as the candidate drug molecules, and 3D structures of their physiologic targets (SERT and AChE) were retrieved from PDB (www.rcsb.org) with designation codes 5I6X (19), 5I71 (19), 5I73 (19), 4M0E (20), 4EY7 (21) and 5HF9 (22).
Molecular docking
Ligandflexible molecular docking simulations were performed with Lamarckian Genetic Algorithm (LGA) (23) incorporated into AutoDock 4.2 (7). All the simulation procedures were conducted according to the previous studies (24). Drugtarget inhibition constants (K_{i}) were estimated through Equations 1 and 2:
(1)
(2)
In Equation 1, ΔG_{b} represents free binding energy (cal.mol^{1}), R is the gas constant (cal.K^{1}.mol^{1}), and T stands for temperature in kelvin (For docking simulation: 298.15 K), and 2.71828 is indicative of a Napier’s constant. In the case of equation 2, E_{vdW}, E_{Hbond}, E_{Desolvation, }and E_{Electrostatic} represent van der Waals energy, hydrogen bond energy, and desolvation energies for drugtarget interaction, and ΔG_{Torsional }is the estimated loss of torsional free energy upon binding to the target.
Experimental design
All statistical analysis and modeling procedures were performed via the BoxBehnken method incorporated into DesignExpert softwarev.7 (StateEase, Corp., and Minnesota) (25). Methodological factors and their assigned levels to construct models are summarized in Table 2. Three levels were considered for each factor under study. Codes were indicative of low (1), medium (0) and upper (+1) levels of the factors, respectively. Appropriate factors and their assigned levels were determined in a way that a broad experimental domain within reasonable endpoints could be scanned.
The subsequent step included the design of experiments (DOE) to offer a BoxBehnken matrix that comprised various docking trials (solutions). Each trial contained different combinations of factor levels. Citalopram and donepezil were docked into the binding sites of SERT and AChE according to DOE trials. A typical matrix for 6 independent factors, each defined in 3 levels, offered 54 docking trials. Results of all docking trials were translated into docking accuracy via Equation 3:
(3)
R_{Theoretical} is the theoretical response or estimated target inhibitory constant (pk_{i,}_{insilico}), whereas R_{Experimental} stands for experimental response or invitro target inhibitory constant (pk_{i,}_{invitro}) driven from standard databases (Binding MOAD, PDB Bind and Binding DB). To explain more, the final endpoint was considered an easily detectable parameter pKi, indicative of drug binding affinity. Secondorder polynomial functions were used within Equation 4 to correlate the responses with designated factors:
(4)
In the above equation, y is the predicted response; β_{0} is an intercept term; βi, βij and βii are linear, quadratic and interaction coefficients, respectively, x_{i} and x_{j }are independent variables in coded levels (1 to 1). The ε value shows a random error. The results were reported by using probability value (pvalue) with 0.05 as the confidence level. Analysis of variance (ANOVA) was implemented for each endpoint to determine the significant factors of the developed model. The models/factors were recognized as significant in each case if the probability value (pvalue) was less than 0.05. Model simplifications were carried out via the elimination of nonsignificant terms (p > 0.05) in all of the model equations. Approved models were characterized by Fvalue, lack of fit Fvalue, predicted Rsquared, adjusted Rsquared and Adeq precision to ensure that they could successfully scan the design space.
Results and Discussion
Internal validation
The validity of the AutoDock4.2 method for docking of selected CNS drugs inside their targets was interpreted in terms of RMSD of ligand atoms in redocked and crystallographic conformations (Table 3). According to the results, all the crystallographic files could pass the filter since AutoDock4.2 could successfully predict the crystallographic (bioactive) conformation (23) within 50 independent GA runs and 2.5 × 106 maximum number evaluations incorporated into LGA. In confirmation of the obtained results, 3D schematic representations of validation results with the best RMSD poses for each drug is depicted in Figure 2.
Model development and statistical analysis
Threelevel BoxBehnken designs are generated by combining twolevel factorial designs and incomplete block designs (26). The technique brings about a few benefits, such as desirable statistical properties and, most importantly, the requirement for only a fraction of trials needed for a 3level factorial. Run number of BoxBehnken design can be estimated according to Equation 5:
(5)
k is the factor number and cp is the replicate number of the central point. In a cubic scheme of BoxBehnken design (Figure 3), a model consists of a central point and the middle points of the edges.
Citalopram
The ANOVA results for matrix responses (ΔpK_{i} 2.2984.111) are summarized in Table 4. Statistical analysis proved the quadratic polynomial model to be highly significant (pvalue < 0.0001) for data fitting. The acquired model in terms of coded values is illustrated by Equation 6:
(6)
As referred, for citalopram, the quadratic model was capable of describing the relation between ΔpK_{i} as a dependent variable and factors B (AutoGrid space), E (Drug optimization method), and F (Target flexibility) as independent variables. In the quadratic model, factors and secondorder interferences with pvalues larger than 0.05 were eliminated by stepwise selection. The lack of fit of the model (1.14) implied that it was not significant with regard to the pure error. Pred Rsquared was estimated to be 0.9789, and moreover, it was in reasonable agreement with Adj Rsquared (0.9842). A good correlation between factors and responses could be confirmed by the Adj Rsquared value, and it meant that most of the variations of response were predictable by model. Adequate precision measures the signaltonoise ratio and the estimated amount (43.939) was indicative of an adequate signal (A ratio greater than 4 is desirable). On the basis of such model characteristics, it was deduced that obtained model could navigate the design space.
According to ANOVA results (summarized in Table 4), ΔpK_{i} sensitivity to the effective factors could be ranked as E>B>>C>D>A>F, and factors C, D, A & F were detected as insignificant model terms (pvalue > 0.1). Factors E (drug optimization method) and B (grid spacing) were significant model terms, while factor E (Fvalue 2204.33) exhibited extremely significant performance. High interactive effects of factors E and F have been observed on the response (pvalue < 0.0001). Quaternion degrees for drug (factor C) was recognized as an insignificant factor in docking of citalopram into SERT. This indicated that docking accuracy was not dependent on the flipping angle of the citalopram molecule. Lack of significant sensitivity toward variations of factor A (torsion degrees for drug) may be interpreted by the fact that docking simulations are commonly initialized by the cocrystallographic conformations, where fitted binding pose of the drug is applied to run the docking procedure. Such a result may have a practical outcome in docking studies; to achieve a desirable result within reasonable computation times, one might set torsions degrees for the ligand at larger values for more rigid structures.
Drug optimization method
Factor E (drug optimization method) was estimated to be the most significant model term in docking of citalopram into ST binding site. The observed highly significant effect might be attributed to the chiral center of citalopram and hence its determinant role in predocking conformation of the drug. At first glance, equation 6 indicated that higher docking accuracies could be expected from AM1based optimization of the citalopram structure method, but due to the highly significant interactive effect of E and F, the best combination toward lowest ΔpK_{i} (2.298) was achieved by PM3 optimization method (coded level of +1). Such a result is in accordance with the inversion barriers of trivalent nitrogen in nitrogencontaining compounds, which are commonly low for AM1 and high for PM3. An apparent consequence is that some nitrogen geometries may be predicted to be flat by AM1 and pyramidal by PM3. Hence it seemed that PM3 could represent a relatively appropriate description of nitrogen geometry and hence more realistic binding interactions with the target.
A common belief is that, unlike ligandbased drug design techniques in which the initial geometry of a bioactive molecule is important, structurebased approaches such as docking are not seriously dependent on the primary optimization of ligand. Indeed beginning a docking practice with a cocrystallographic binding pose of a ligand is a common approach. A rationale is that during molecular docking simulations, molecular conformations are varied via changes in torsion, translation, and quaternion. But the different scenario that was observed with the present study was the highly significant effect of the optimization method (Factor E) on docking output. On the basis of obtained results, it may be assumed that chiral molecules, particularly those bearing nitrogen atoms within a nonpolar scaffold, can undergo an appropriate semiempirical method such as PM3 to afford better results.
Target flexibility
Target flexibility was incorporated into our modeling study via considering different holostructures of SERT. The results of the statistical analysis were in accordance with what we expected. Higher docking accuracies could be attained via docking of citalopram into the binding site of the SERT structure with the highest resolution (PDB code 5I6X). It was found that decreasing the resolution of SERT from 3.14 Å (5I6X) to 3.24 Å (5I73) reduced docking accuracies (Figure 9). However,, future studies may be directed toward selecting more induced fit models of the protein and statistical analysis through central composite design (CCD techniques).
Grid spacing
A grid map comprises a 3D frame of regularly spaced points for incorporating the target. On the basis of ligand atom types, a probe atom corresponding to the atom type is placed at each grid point, and the energy of interaction of each probe atom (grid point) with surrounding macromolecular atoms is estimated and assigned to the corresponding grid point (24) (Figure 5). Grid spacing (Factor B) is designated as the distance between adjacent grid points. Grid spacing is the distance between adjacent AutoGrid points. ANOVA showed that if the grid spacing is set to lower values, higher AutoDock accuracies for the citalopramSERT complex will be achieved in confirmation of previous results (13). Lower grid spacings increase the precision of probe scanning within the designated grid box, and this would probably be translated into better SERT inhibition constants.
Interactive factor effects
Factor interaction is likely to occur whenever different responses are generated on the basis of different settings of two factors. This dependence of factor levels to each other may be best interpreted by interaction plots (Figure 6). In this case, interactive factors will be depicted by two nonparallel lines, implying that the effect of one factor depends on the level of the other. ANOVA results proved highly significant interactive effects between factors E (drug optimization method) and F (target flexibility) with pvalue < 0.0001.
It was indicated that the significant effect of factor E on estimated SERT inhibition constants was more pronounced at lower levels of F. As could be seen from the graph in Figure 6, the red line is indicative of the effect of the drug optimization method (F) within a SERT 3D structure with PDB code 5I6X and the black line represents the effect of drug optimization method (F) within a SERT 3D structure with PDB code 5I73. Higher docking accuracies might be expected when other factors (A, B, C & D) were held at their lower levels (such as factor F) (Figure 6).
Interaction plots displayed a cross point, and the location of this point showed a distinctive situation within model space in which relatively similar SERT inhibition constants could be expected by docking into all PDB driven 3D SER structures (levels 1 & +1) if the cocrystallographic conformation of citalopram is set as the starting point (Midlevel of factor E).
The 3D surface is known as the “response surface” provided a perspective visualization of factor effects on the response at different levels of other factors. 3D response plots were developed to indicate the simultaneous effect of interactive term EF on docking accuracy (Figure 7). In confirmation of our previous results, a surface was steep and indicated that the interaction between two factors was highly significant. More accurate SERT inhibition constants might be predicted by running the PM3 semiempirical method (higher levels of factor E) at declined levels of other factors.
Numerical optimization
DOE provides a series of solutions (optimum combinations of factor levels) to achieve the most desirable responses. For this purpose, optimization criteria for levels of factors A, B, C, and D were set in range (spanning from 1 to +1) while factors E (drug optimization method) and F (target flexibility) were set at precise levels 1, 0 and +1 since they were categorical but not numerical factors. With the aim of achieving optimized solutions, the goal for response (ΔpK_{i}) was primarily set at a minimum.
It should be noted that in each case, solutions with desirability equal to 1 were picked up as optimum. Desirability is an objective function ranging from 0 (worst condition) to 1 (ideal case). This function transforms each response value to a desirability index. The program looks for the largest desirability index and presents a series of solutions that best maximize the desirability index. Obtained results showed that the most accurate predictions of SERT inhibition constant (minimum ΔpK_{i}) could be envisaged through various simulations conditions, and careful selection of factors led to highly enhanced accurate responses (Table 5). However, it should be emphasized that choosing the best solution depends on financial and time restrictions. A characteristic feature in all of the proposed docking solutions is the lower level of factor F and a higher level of factor E, which confirmed the previous results of this study.
Donepezil
In the case of donepezil, ANOVA results for the responses (ΔpK_{i} 0.2405.465) are summarized in Table 6. Statistical analysis proved the quadratic polynomial model to be highly significant with a low probability value (pvalue < 0.0001) for data fitting (Equation 7 in terms of coded values).
(7)
Lack of fit Fvalue (1.74) implied that lack of fit was not significant with regard to the pure error. There is a 13.36% chance that a “Lack of Fit Fvalue” this large could occur due to noise. Pred Rsquared (0.9820) was in reasonable agreement with Adj Rsquared (0.9888). A good correlation between factors and responses could be confirmed by the Adj Rsquared value. Adequate precision (43.939) was indicative of an adequate signal. On the basis of obtained data, it was deduced that model could navigate the design space.
According to ANOVA results (summarized in Table 6), factor effects could be ranked as B>E>>D>C>A>F while D, C, A & F were insignificant model terms (pvalue > 0.1). It was found that factors B (grid spacing) and E (drug optimization method) were significant model terms. Among the pairwise interactions, EF was the significant model term (pvalue<0.0001) followed by BE (pvalue 0.0102).
Grid spacing
Factor B (grid spacing) was the most significant model term. The polynomial quadratic model (Equation 7) predicted better AChE inhibition constants for donepezil at shorter grid spacings (0.3 Å). The effect was similar to citalopram but more noticeable in the case of donepezil. One possible explanation is the presence of bulky molecular structure of donepezil that necessitates shorter grid spacings in docking simulations.
Drug optimization method
Algebraic signs of quadratic model terms (Equation 7) indicated that the application could expect higher docking accuracies of the PM3 method for primary optimization of donepezil structure. Comparative statistical inspection of the results showed that with regard to pvalues, although being significant, the effect of factor E is more significant for citalopram. This may be attributed to the following explanations:
Unlike citalopram, donepezil includes one nitrogen atom within a hydrophobic structural pattern and inversion barriers of trivalent nitrogen for AM1 and PM3 semiempirical methods, less dependent on the PM3 optimization method might be explainable.
More flexible structure of citalopram (more active torsions) with regard to donepezil.
Interactive factor effects
ANOVA results demonstrated a significant pairwise interactive effect between factors E and F (pvalue < 0.0001). The significant effect of factor E was most pronounced when donepezil was docked into the AChE model that possessed the highest resolution (PDB 4M0E) (Figure 8). The observed interaction pattern was different from that of citalopram. More accurate enzyme inhibitory activities for donepezil could be expected within two scenarios; AM1based optimization of drug molecules and docking simulations on 4EY7 or PM3based optimization of a drug molecule with docking simulations on 4M0E.
All the interaction plots of EF interactive effects displayed a cross point on midlevels of factor E (Initial Cocrystallographic conformation of donepezil). To explain more, when cocrystallographic conformation of donepezil was used as the starting point for docking simulations, the estimated AChE inhibition constant was not seriously dependent on the selected PDB model of the target. Such interferences might not be detected via applying one factor at each time method.
Onefactor plots confirmed the direction of interactive effects and indicated the highly significant effect of factor E was detected when other factors were held at their upper levels.
3D plots representing simultaneous effects of factors B and E at different levels of other factors are depicted in Figure 10. In upper levels of other factors, more desirable docking accuracies were expected at higher levels of E (PM3 or PM3like optimization methods). The surface in midlevels of factors is relatively smooth that was indicative of a less significant interactive effect between B and E.
3D response plots were also developed to interpret the interactive EF effect on docking accuracy. As could be seen from the plots (Figure 11), docking accuracy tended to increase at higher levels of factor E as the levels of other factors declined to lower levels. 3D plots obviously showed that when factor levels were held at their midlevels, no desirable docking accuracy would be expected. 3D surface plots in lower levels of factors A, C, D and F showed that desirable docking accuracies could be attained at lower levels of factor E and any level of factor B.
3D response plots were also applied to indicate the simultaneous effect of interactive term EF on docking accuracy (Figure 11). It was found that responses to factor levels fitted a hyperbolic pattern with relative symmetric distribution and steep surfaces. This could be related to the highly significant effect of EF on the response (pvalue < 0.0001), which is not seriously dependent on the levels of other factors. As could be seen from the plots, more reliable results may be assumed at lower levels of both E and F or higher levels of both E and F. Such interaction pattern can be demonstrated that when higher resolution PDB conformation of AChE is used for docking of donepezil, it would be better to optimize the drug structure with PM3 method while the reverse is true when the lower resolution of AChE conformation is applied. This was also previously confirmed by the interaction plots.
Numerical optimization
All the optimization criteria for factors A, B, C, D, E, and F were set as before, and the goal for docking accuracy (ΔpK_{i}) was fixed at a minimum. On the basis of offered optimized solutions, maximum docking accuracy (minimum ΔpK_{i}) might be achievable via various conditions (Table 7). However, choosing the best solution depends on the financial and time limitations.
In confirmation of ANOVA results, it was revealed that the most accurate predictions of AChE inhibition constant (minimum ΔpK_{i}) could be envisaged when both of the factors F and E were set at their upper or lower levels.
Characteristics of candidate CNS drugs and their mechanisms of pharmacological action
Actual/coded values of selected factors for AutoDock4.2 based RSM study of citalopram and donepezil
Factors understudy  Factor levels  

Low: Actual (Coded)  Medium: Actual (Coded)  High: Actual (Coded)  
A: Torsion degrees for drug  5 (1)  20 (0)  50 (+1) 
B: Grid spacing (Å)  0.3 (1)  0.375 (0)  0.5 (+1) 
C: Quaternion degrees for drug  5 (1)  20 (0)  50 (+1) 
D: Translation (Å)  0.2 (1)  0.3 (0)  0.5 (+1) 
E: Drug optimization method  AM1 (1)  Cognate ligand (0)  PM3 (+1) 
F: Target flexibility  Lowest resolution (Å)  Medium resolution (Å) PDB code (0)  Highest resolution (Å) PDB code (+1) 
AutoDock 4.2 validation results for different holo PDB structures of intended CNS targets
Drug/No.  PDB ID  Resolution of  Number of GA runs  No. of conformationsin topranked cluster  Maximum  RMSD from Reference (Å) 

Citalopram  
1  5I6X  3.14  50  49  2.5 × 10^{6}  1.58 
2  5I71  3.15  50  49  2.5 × 10^{6}  0.90 
3  5I73  3.24  50  50  2.5 × 10^{6}  0.82 
Donepezil  
1  4M0E  2.00  50  50  2.5 × 10^{6}  0.55 
2  5HF9  2.20  50  50  2.5 × 10^{6}  1.92 
3  4EY7  2.35  50  47  2.5 × 10^{6}  0.46 
ANOVA report for significant terms of a polynomial quadratic model in docking study of citalopramserotonin transporter complex
Source  Sum of squares  df  Mean square  Fvalue  Pvalue 

Model  19.41  11  1.76  301.22  <0.0001 
B  0.082  1  0.082  14.02  0.0005 
E  12.91  1  12.91  2204.33  <0.0001 
EF  5.59  1  5.59  954.33  <0.0001 
Residual  0.25  42  5.858E003  
Lack of fit  0.15  25  6.160E003  1.14  0.3987 
Pure error  0.092  17  5.415E003  
Cor Total  19.66  53 
RSMbased optimum solutions for AutoDock4.2 simulations (in terms of coded factor levels) leading to the most accurate inhibition constants of serotonin transporter by citalopram (1 > ΔpKi); A: Torsion degrees for drug; B: Grid spacing (Å); C: Quaternion degrees for drug; D: Translation (Å); E: Drug optimization method; F: Target flexibility
No. optimized solution  Factor levels  

A  B  C  D  E  F  ΔpK_{i}  
1  0.33  0.96  0.85  0.14  1.00  1.00  0.925 
2  0.42  0.79  0.67  0.50  1.00  1.00  0.933 
3  0.25  0.73  0.57  0.85  1.00  1.00  0.934 
4  0.42  0.40  0.90  0.85  1.00  1.00  0.938 
5  0.97  0.41  0.78  0.04  1.00  1.00  0.960 
6  0.85  0.93  0.67  0.27  1.00  1.00  0.968 
7  0.69  0.94  0.36  0.84  1.00  1.00  0.987 
8  0.32  0.40  0.04  0.73  1.00  1.00  0.996 
9  0.70  0.50  0.28  0.32  1.00  1.00  0.999 
ANOVA results for significant terms of a polynomial quadratic model in docking study of donepezilAChE complex
Source  Sum of squares  df  Mean square  Fvalue  Pvalue 

Model  248.19  20  12.41  234.38  <0.0001 
B  1.19  1  1.19  22.42  <0.0001 
E  0.91  1  0.91  17.18  0.0002 
BE  0.39  1  0.39  7.42  0.0102 
EF  189.12  1  189.12  3571.88  <0.0001 
Residual  1.75  33  0.053  
Lack of fit  1.08  16  0.068  1.74  0.1336 
Pure error  0.66  17  0.039  
Cor Total  249.94  53 
RSMbased optimum solutions for AutoDock4.2 simulations (in terms of coded factor levels) leading to the most accurate inhibition constants of AChE by citalopram (0.2 > ΔpKi); A: Torsion degrees for drug; B: Grid spacing (Å); C: Quaternion degrees for drug; D: Translation (Å); E: Drug optimization method; F: Target flexibility
No. Optimized solution  Factor levels  

A  B  C  D  E  F  ΔpK_{i}  
1  0.45  0.99  0.93  0.97  1.00  1.00  0.152 
2  0.99  1.00  0.90  0.99  1.00  1.00  0.160 
3  0.06  0.86  1.00  0.91  1.00  1.00  0.161 
4  0.08  0.49  1.00  1.00  1.00  1.00  0.187 
5  0.11  0.89  0.98  0.96  1.00  1.00  0.185 
Conclusion
Availability of facile, timeefficient, and accurate computeraided or insilico drug design techniques is an urgent requirement for identifying and developing potent and selective medicinal agents. Within the structurebased strategies, molecular docking is a frequently used and valuable computational method for matching ligands/drugs into the environment of a validated target. Despite several advantages and fruitful historical outcomes, current docking simulations are mostly restricted to inaccurate estimated binding affinities. In light of the above explanations, improving docking accuracy to fill the gap between theoretical and experimental data through statistical optimization of effective variables may be plausible. Efficient statistical techniques such as RSM can be appropriately utilized for the identification of effective factors and their optimization toward more robust docking simulations. RSMs offers a substantial advantage over commonly applied onefactorateachtime techniques in a way that, besides individual factor effects, interactive effects may also be considered within noticeably fewer trials. Within the present contribution, the full potential of RSM in optimizing molecular docking simulations was unveiled through BoxBehnken derived ANOVA analysis of AutoDock4.2 based binding affinity prediction. For this purpose, polynomial, quadratic models were constructed for the binding of highly prescribed antidepressant (citalopram) (R^{2} 0.9789) and antiAlzheimer’s (donepezil) (R^{2} 0.9820) drugs to physiological targets SERT and AChE. Significant individual and interactive factor effects on the accuracy of estimated target inhibition constants were statistically elucidated. It was revealed that estimated binding affinities citalopram and donepezil were mostly affected by the predocking optimization method and AutoGrid spacing, respectively. One of the advantageous features of RSMs is the identification of interactive effects that simultaneously change the response. For citalopram, the optimization method exhibited significant pairwise interaction with conformational flexibility of SERT, while in the case of donepezil, the binding of the drug to AChE was significantly affected by the interactive effect of grid spacing with the optimization method. Probably the most productive section of study results was the numerical optimization that offered a few optimized docking simulations leading to significantly higher accuracies in AutoDock4.2 driven SERT and AChE inhibition constants. The outputs of this study may indicate the full potential of RSMs for the development of optimized AutoDock protocols toward the rational design of privileged medicinal scaffolds.
Acknowledgements
References

1.
Domínguez A, Álvarez A, Hilario E, SuarezMerino B, GoñideCerio F. Central nervous system diseases and the role of the bloodbrain barrier in their treatment. Neurosci. Discov. 2013;1:314.

2.
Wittchen HU, Jacobi F, Rehm J, Gustavsson A, Svensson M, Jonsson B, Olesen J, Allgulander C, Alonso J, Faravelli C, Fratiglioni L, Jennum P, Lieb R, Maercker A, van Os J, Preisig M, SalvadorCarulla L, Simon R, Steinhausen HC. The size and burden of mental disorders and other disorders of the brain in Europe 2010. Eur. Neuropsychopharm. 2011;21:65579.

3.
Jindal H, Bhatt B, Sk S, Malik JS. Alzheimer disease immunotherapeutics: Then and now. Hum. Vacc. Immunother. 2014;10:27413.

4.
Mandal S, Moudgil M, Mandal SK. Rational drug design. Eur. J. Pharmacol. 2009;625:90100. [PubMed ID: 19835861].

5.
Dar AM, Mir S. Molecular Docking: Approaches, types, applications and basic challenges. J. Anal. Bioanal. Tech. 2017;8:3569.

6.
Kuntz ID, Blaney JM, Oatley SJ, Langridge R, Ferrin TE. A geometric approach to macromoleculeligand interactions. J. Mol. Biol. 1982;161:26988. [PubMed ID: 7154081].

7.
Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS, Olson AJ. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comput. Chem. 2009;30:278591. [PubMed ID: 19399780].

8.
Jones G, Willett P, Glen RC, Leach AR, Taylor R. Development and validation of a genetic algorithm for flexible docking. J. Mol. Biol. 1997;267:72748. [PubMed ID: 9126849].

9.
Rarey M, Kramer B, Lengauer T, Klebe G. A fast flexible docking method using an incremental construction algorithm. J. Mol. Biol. 1996;261:47089. [PubMed ID: 8780787].

10.
Friesner RA, Banks JL, Murphy RB, Halgren TA, Klicic JJ, Mainz DT, Repasky MP, Knoll EH, Shelley M, Perry JK, Shaw DE, Francis P, Shenkin PS. Glide: A new approach for rapid, accurate docking and scoring Method and assessment of docking accuracy. J. Med. Chem. 2004;47:173949. [PubMed ID: 15027865].

11.
Khuri AI, Mukhopadhyay S. Response surface methodology. WIREs Comp. Stats. 2010;2:12849.

12.
Ghasemi N, Bandehpour M, Ranjbar J. Optimization of key factors in serum free medium for production of human recombinant GMCSF using response surface methodology. Iran. J. Pharm. Sci. 2019;18:14656.

13.
Bohlooli F, Sepehri S, RazzaghiAsl N. Response surface methodology in drug design: A case study on docking analysis of a potent antifungal Fluconazole. Comput. Biol. Chem. 2017;67:15873. [PubMed ID: 28129567].

14.
Kracmarová A, Drtinová L, Pohanka M. Possibility of acetylcholinesterase overexpression in Alzheimer disease patients after therapy with acetylcholinesterase inhibitors. Acta Med. 2015;58:3742.

15.
Apostolova LG. Alzheimer’s disease. Continuum. 2016;22:41934. [PubMed ID: 27042902].

16.
Sahraian S, Babashams M, RezaSoltani P, Najmabadi H, Kahrizi K, Gorgani SH. Serotonin transporter polymorphism (5HTTLPR) and citalopram effectiveness in Iranian patients with major depressive disorder. Iran. J. Psychiatry. 2013;8:8691. [PubMed ID: 24130607].

17.
Cipriani A, Purgato M, Furukawa TA, Trespidi C, Imperadore G, Signoretti A, Churchill R, Watanabe N and Barbui C. Citalopram versus other antidepressive agents for depression. Cochrane Database Syst. Rev. 2012;7:CD006534.

18.
AgatonovicKustrin S, Kettle C, Morton DW. A molecular approach in drug development for Alzheimer’s disease. Biomed. Pharmacother. 2018;106:55365. [PubMed ID: 29990843].

19.
Coleman JA, Green EM, Gouaux E. Xray structures and mechanism of the human serotonin transporter. Nature. 2016;532:3349. [PubMed ID: 27049939].

20.
Cheung J, Gary EN, Shiomi K, Rosenberry TL. Structures of human acetylcholinesterase bound to dihydrotanshinone I and territrem B show peripheral site flexibility. ACS Med. Chem. Lett. 2013;4:10916. [PubMed ID: 24900610].

21.
Cheung J, Rudolph MJ, Burshteyn F, Cassidy MS, Gary EN, Love J, Franklin MC, Height JJ. Structures of human acetylcholinesterase in complex with pharmacologically important ligands. J. Med. Chem. 2012;55:102826. [PubMed ID: 23035744].

22.
Franklin MC, Rudolph MJ, Ginter C, Cassidy MS, Cheung J. Structures of paraoxoninhibited human acetylcholinesterase reveal perturbations of the acyl loop and the dimer interface. Proteins. 2016;84:124656. [PubMed ID: 27191504].

23.
GhorbanDadrass L, Madadkar SA, Shafiei A, Mahmoudian M. Flexible ligand docking studies of matrix metalloproteinase inhibitors using lamarkian genetic algorithm. DARU. 2004;12:110.

24.
RazzaghiAsl N, Ebadi A, Edraki N, Shahabipour S, Miri R. Ab initio modeling of a potent isophthalamidebased BACE1 inhibitor: Amino acid decomposition analysis. Med. Chem. Res. 2013;22:325969.

25.
Box GEP, Draper NR. Empirical modelbuilding and response surfaces. 87th Ed. New York: John Wiley & Sons, Inc; 1987.

26.
Ferreira SL, Burns RE, Ferreira HS, Matos GD, David JM, Brandão GC, da Silva EG, Portugal LA, dos Reis PS, Souza AS, dos Santos WN. BoxBehnken design: An alternative for the optimization of analytical methods. Anal. Chim. Acta. 2007;597:17986. [PubMed ID: 17683728].