The life cycle of the virus as many other viruses has an essential proteolytic auto-processing step. 2 3-4 The 3CL-protease takes an important role in freeing individual functional proteins from single-chain polyprotein that has been translated by the host cell translation machinery. 5 Interfering to this step has been shown to have capability of inhibiting virus replication. 6 Therefore, 3CL-protease is one of the valid and most attractive drug targets for Covid-19
treatment.
Developing a drug molecule from scratch for curing a disease is quite expensive and time consuming. Instead, repurposing an FDA approved drug is a much faster and less expensive strategy since the time is very critical parameter in controlling such invasive viruses. Even testing all FDA approved drugs experimentally one by one may cause intolerable waste of the time. To help to ease the problem, computational calculations can be very handy. There have already been several computational studies about predicting an efficient FDA approved drugs that can be repurposed for COVID-19. [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] In searching for drug candidates as inhibitors, there are two core questions that can be addressed: (a) what is the binding mode? (i.e. where does it bind from the enzyme?) (b) what is the binding affinity? (binding free energy). 18 Inhibition of an enzyme can be involved in various mechanisms. For instance, an inhibitor can bind to the enzyme from the active site in the absence of substrate (competitive inhibition) forming an additional or it can bind in the presence of the substrate (i.e. enzyme-substrate complexed) from the allosteric site (uncompetitive inhibition).
Alternatively, it can reduce the activity of the enzyme by binding from allosteric site of the enzyme either free or complex (non-competitive inhibition). 19 In all these mechanisms the equilibrium constant, K s belonging to the reaction between the enzyme and substrate is shifted towards the reactants by existing inhibitor in the media. Kinetic studies to understand the activity loss of the enzyme on substrate by the inhibitor report values such as inhibition equilibrium constant, K i or half maximal inhibitory concentration, IC50. Both these experimental parameters are proportional (they are equal in special cases). Thermodynamic equilibrium constant K i is related to the standard binding Gibbs free energy of the system by:
The ΔG o can also be calculated from the thermodynamic potentials. Thus, it is possible to compare experimental ΔG o determined from K i and theoretical ΔG o predicted by thermodynamic potentials using computational methods. However, logarithmic relation between the K i and ΔG o complicates the reliability of the calculations. An error of 3-4 kcal is acceptable for the most computations, but this makes thousands of times less efficient inhibitor. Therefore, most computational methods fail to estimate the experimental binding affinity trend (relative binding free energies) among the inhibitors. Reproducing the absolute standard free energies of binding is even more difficult for the computational methods. Moreover, a trade-off between the computational accuracy and speed must be made.
However, almost all of the studies screen databases listing the FDA approved drugs at the molecular docking level, which is mostly based on geometrical alignment of a ligand into a binding pocket, namely the active site of 3C-like protease (3CL pro ). Despite the speed of docking approach, they are very coarse methods and almost never predict the correct experimental binding affinity trend among the inhibitors. [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] More sophisticated methods to calculate the potential binding free energy of inhibitor candidate to the protein ranges from post molecular dynamics simulations such as Molecular
Mechanics Poisson-Boltzmann Surface Area (MMPBSA) [20] [21] [22] [23] to perturbation methods such as Bennett acceptance ratio (BAR), 24-31 the latter being much more accurate. [42] [43] [44] are also commonly used and proven to be sufficiently accurate in estimating solvation free energies along with relative and absolute binding free energies. However, all these methods require very high computational cost. Screening all the drug candidates using these methods may not be very practical. For a typical free energy of binding, these methods require at least 30 λ-states for the decoupling of L from PL complex and another that many for decoupling L from water (i.e. solvation free energy of ligand), which increases the computational cost by 60-folds.
Machine Learning (ML) techniques have attracted a great interest for the past decade due to their successful algorithms to scientific questions including but not limited to chemical reactions, [45] [46] potential energy surfaces, [47] [48] [49] [50] forces, 51-53 atomization energies, 54-55 and proteinligand complex scorings. 56 . One of the most promising aspects of the ML techniques is that the trained models can be applied to new systems (transferrable). 57 Recently several models using active learning such as ANI-1, 58 ANI-1x 59 and ANI-1cxx 60 Here, we introduce a new strategy to estimate free energies of solvation of small organic compounds and binding to proteins in explicit solvent using single end-state MD simulations.
The method is adopted from ANI-1ccx neural network potentials (Machine Learning) for the The insertion of the ligand to an environment of solvent (solvation free energy) or receptor (binding free energy) can be defined by a coupling parameter, λ. 61 At each λ-state, the potential energy of the system would be:
where V is the interaction energy between the ligand and environment.
ܷ is the reference potential energy, which includes all energy terms except for the ligand-environment interaction.
The partition function is;
where F is the free energy depending on λ-state. Re-organizing the above equation,
The equation above is only valid when the ligand-environment interaction is proportional to coupling term λ. This results in;
In the limiting case, this equation yields as;
The equation implies that the free energy change is nearly half-the potential between the ligand-environment interaction in the complex configurational ensemble. 20 For the binding studies, the protein-ligand complex (~4700 atoms) was placed in the center of a dodecahedron box. For the solvation studies, the ligand was placed in the center of a cubic box.
Each system was solvated in the TIP3P model type water 64 Energy minimization was carried out to a maximum 100 kJ.mol -1 .nm -1 force using Verlet cutoff scheme. For both long range electrostatic and Van der Waals interactions, a cutoff length of 12 Å was used with the Particle Mesh Ewald method (PME) (6 th order interpolation). 65 The neighbor list update frequency was set to 1 ps -1 . As with our earlier studies, [66] [67] a stepwise energy minimization and equilibration schemes were used. Each minimization step consisted of a 5000 cycle of Steepest Descent and a subsequent 5000 cycles of l-bfgs integrators.
After minimization, each system was equilibrated within three steps using Langevin In order to test the ML-QM method, we run 10-ns MD simulations for the molecules whose solvation free energies have been extensively studied. We first built a system in which the ligand is replaced in the center of box and the box is solvated with just single water molecule. The ligands studied are given in Figure 1 . We prepared other independent simulations in which we increased the number water molecules (nH 2 O) in the box. Each time the average ML-QM energies were almost 0. This has been the case until nH 2 O=100. Figure 3a shows the interaction energy between the ligand and water molecules calculated by ML-QM energy difference:
is the total energy of the system including ligand and water.
are the energies of free ligand and free water molecules, respectively, which have been calculated by extracting from the trajectory of complex structure. When we kept adding more water in the system and running other MD simulations, we found that starting from nH 2 O=100, the system undergoes from the unbound state to bound state by a drastic energy jump after ~9.5 ns. The degree of freedom for the water molecules in the bound state must be It is apparent that the ML-QM interaction energy is constant for the simulations with more than sufficient waters exist. The empirically determined linearity constant is α =1/2, thus the free energy change is approximately half of the interaction energy between the ligand and environment. However, considering the experimental conditions are more or less similar, the energy profile (i.e. the trend among the inhibitors) should at least be reproduced. Here, we see that ML-QM is the most accurate method to yield the correct trend among the compared methods ( Figure 5 ). 