Keywords: de novo design, drug discovery, protein inhibitors, ligand growth, combinatorial chemistry, HIV-1 protease

We present a new open-source software, called OpenGrowth, which aims to create de novo ligands by connecting small organic fragments in the active site of proteins. Molecule growth is biased to produce structures that statistically resemble drugs in an input training database. Consequently, the produced molecules have superior synthetic accessibility and pharmacokinetic properties than randomly grown molecules. The growth process can take into account the flexibility of the target protein and can be started from a seed to mimic R-group strategy or fragment-based drug discovery. Primary applications of the software on the HIV-1 protease allowed us to quickly identify new inhibitors with a predicted Kd as low as 18nM. We also present a graphical user interface that allows a user to select easily the fragments to include in the growth process.

Introduction

Computational tools are now fully integrated in all steps of structure-based drug discovery, from the early stage of homology modelling [START_REF] Cavasotto | Homology modeling in drug discovery: current trends and applications[END_REF] to free energy perturbations used during lead optimization [START_REF] Mobley | Perspective: Alchemical free energy calculations for drug discovery[END_REF] or predictions of pharmacokinetic properties [START_REF] Norinder | Prediction of ADMET Properties[END_REF] . High-throughput virtual screenings (and especially molecular docking) have proven to be reliable tools to identify new hits [START_REF] Kitchen | Docking and scoring in virtual screening for drug discovery: methods and applications[END_REF] and are widely used in drug discovery.

Once a target protein has been identified, this kind of approach requires a predefined library of compounds to screen. The most straightforward libraries consist of existing molecules, either commercially available ones (such as the ZINC library [START_REF] Irwin | ZINC--a free database of commercially available compounds for virtual screening[END_REF] ), a set of natural products, in-house compounds or current drugs/drug candidates. This methodology leads either directly to new candidates or to hits that will later be optimized. The former, a "drug repurposing" approach, has proven to be both efficient and effective [START_REF] Mullard | Drug repurposing programmes get lift off[END_REF] . A prime example is raloxifen, which is now used to treat both osteoporosis and breast cancer [START_REF] Riggs | Selective estrogen-receptor modulators --mechanisms of action and application to clinical practice[END_REF] .

However, drug repurposing is limited in scope, especially when we aim to find original scaffolds.

Increases in computing power and in automation of biochemical assays forecast that large libraries (if not the full ZINC library) will in the near future have been docked or tested in almost all target proteins with known structures. Since the needs for new ways to cure diseases will never stop -because of either drug-resistance [START_REF] Clavel | HIV drug resistance[END_REF] or the emergence of new diseases-, we cannot rely only on already existing molecules or scaffolds in drug discovery. Thereby, there is a need to discover new chemical materials.

So-called computational de novo drug design approaches have emerged to overcome the limitation of existing libraries. In such strategies, a computer generates novel ligands for a target protein.

One must note that in addition to protein ligands, de novo design has also been applied to solar cells [START_REF] Hachmann | The Harvard Clean Energy Project: Large-Scale Computational Screening and Design of Organic Photovoltaics on the World Community Grid[END_REF] and polymers [START_REF] Tew | De Novo Design of Antimicrobial Polymers, Foldamers, and Small Molecules: From Discovery to Practical Applications[END_REF] . Due to the importance of the field, de novo drug design methods have already been reviewed on several occasions [START_REF] Hartenfeller | Enabling future drug discovery by de novo design[END_REF] and a list of available programs and successes of these methods are available in these reviews. We want to present here a brief overview on some of the possible methods to construct the new molecules. First, a common approach is to link two fragments that have been previously docked in the active site (LUDI [START_REF] Bohm | The computer program LUDI: a new method for the de novo design of enzyme inhibitors[END_REF] and PhDD [START_REF] Huang | PhDD: A new pharmacophore-based de novo design method of drug-like molecules combined with assessment of synthetic accessibility[END_REF] are programs that do this). The BREED program recombines molecules by exchanging fragments after the identification of overlapping bonds [START_REF] Pierce | BREED: Generating Novel Inhibitors through Hybridization of Known Ligands. Application to CDK2, P38, and HIV Protease[END_REF] . Another approach, used by FLUX [START_REF] Fechner | Flux (1): A Virtual Synthesis Scheme for Fragment-Based de Novo Design[END_REF] , mutates fragments from bound ligands using retrosynthetic rules such as RECAP [START_REF] Lewell | RECAP--retrosynthetic combinatorial analysis procedure: a powerful new technique for identifying privileged molecular fragments with useful applications in combinatorial chemistry[END_REF] . A fourth strategy -inspired by biology-uses evolutionary algorithms: each moiety of a known ligand is considered a gene that can evolve [START_REF] Lameijer | The molecule evoluator. An interactive evolutionary algorithm for the design of drug-like molecules[END_REF] . Many recently published programs use this method, such as AutoGrow [START_REF] Durrant | AutoGrow: A Novel Algorithm for Protein Inhibitor Design[END_REF] . Starting from a seed, a R-group strategy can also be used: fragments are added around the scaffold to decorate it [START_REF] Pearce | An Automated Workflow for Efficient Structure-Based Lead Optimization[END_REF] . This strategy is not unlike the strategy used by medicinal chemists after the identification of a hit. Finally, growing methods have been successfully implemented, where new fragments (or atoms) are added at a time until a given threshold is reached. This approach shares several similarities to the R-group strategy, as the latter is a growing method with an evolving seed.

Recombination, mutation, evolution and R-group approaches need first to have information about existing ligands, whereas only linking and growing approaches can create completely new ligands. Two sorts of programs using a growing strategy exist. First, it is possible to grow libraries without considering a protein: the produced molecules are then docked by regular means in the target protein.

Thus, there is no need to use three-dimensional representation of molecules: for instance, we described in 2009 an algorithm (called FOG) that analyses a training drug database in terms of probabilities that two fragments are connected [START_REF] Kutchukian | FOG: Fragment Optimized Growth algorithm for the de novo generation of molecules occupying druglike chemical space[END_REF] . New libraries of molecules (that occupy a desired chemical space) are then grown using these probabilities and written as SMILES strings. The second family of program, by far the largest, uses three-dimensional description of the ligands in the active site of a protein. The early approaches were mainly atom-based, i.e. new structures were created atom by atom (e.g. LEGEND [START_REF] Nishibata | Automatic creation of drug candidate structures based on receptor structure. Starting point for artificial lead generation[END_REF] (1991), GenStar 22 (1993) or RASSE 23 (1996)). However, the last notable purely atom-based program was published in 1996 11a, [START_REF] Luo | RASSE: a new method for structure-based drug design[END_REF] and most of the recent algorithms are now fragment-based. In 1993, one of the first fragment-based growth algorithms, GroupBuild [START_REF] Rotstein | GroupBuild: a fragment-based method for de novo drug design[END_REF] , was published. Using a starting fragment (or seed), GroupBuild generates a set of new ligands and randomly keeps one of the top-scored ligands before proceeding to the next step. GrowMol 25 (1994) is both atom and fragment-based and was the first program to use a Metropolis criterion to accept or reject the new fragment. Our group developed the SMoG algorithm [START_REF] Dewitte | SMoG: de Novo Design Method Based on Simple, Fast, and Accurate Free Energy Estimates. 1. Methodology and Supporting Evidence[END_REF] in 1996 whose effectiveness was demonstrated in 2002 [START_REF] Grzybowski | Combinatorial computational method gives new picomolar ligands for a known enzyme[END_REF] . However, the ligands grown by SMoG often lacked "chemical realism" and SMoG was resultantly implemented less than expected. In 2006, BOMB [START_REF] Barreiro | From docking false-positive to active anti-HIV agent[END_REF] and FlexNovo [START_REF] Degen | FlexNovo: Structure-Based Searching in Large Fragment Spaces[END_REF] were released. The details of the BOMB algorithms have never been published and it is thus difficult to compare it with other approaches; FlexNovo was commercial and is no more available for download. LigBuilder was first published in 2000, with a second version presented in 2011 [START_REF] Wang | LigBuilder: A Multi-Purpose Program for Structure-Based Drug Design[END_REF] ; it is described below. In 2012, Vitae Pharmaceuticals described Contour, an internal project that is based on SMoG [START_REF] Ishchenko | Structure-based design technology contour and its application to the design of renin inhibitors[END_REF] ; since it is proprietary software, it is not available to the public.

Thereby, to the best of our knowledge, only one program which uses a growing approach and which is non-proprietary or non-commercial is currently available to the public: LigBuilder (which is free for academic users, but not open-source). This program grows molecules in "growing sites" that must have been pre-defined. Thus, only fragments that fit the shape of the cavity (according to some predefined rules) can be added. In order to increase synthetic accessibility, LigBuilder uses retrosynthetic analysis in combination with connection rules (e.g., O-O bonds cannot be created). We present here a new program that uses other strategies as the ones in LigBuilder. First, we do not define in advance sites of growth: should the placement of a new fragment in a site be not favorable, our approach will reject the fragment (because of energetic criterion or detection of steric clashes). The reason of the existence of the cavities in LigBuilder is to avoid overgrowth: since our program is very fast, we can relax this constraint and our approach does not need a prior preparation of the target protein. Concerning the synthetic accessibility, the use of retrosynthetic analysis or a set of chemical rules to create new molecules is efficient, but can also be too restrictive and may not represent the diversity offered by organic synthesis. For example, the flexibility offered by protection groups is difficult to implement in computational retrosynthetic analysis. Moreover, when a set of chemical reactions is used (such as in SYNOPSIS [START_REF] Vinkers | SYNOPSIS: SYNthesize and OPtimize System in Silico[END_REF] or DOGS [START_REF] Hartenfeller | DOGS: reaction-driven de novo design of bioactive compounds[END_REF] ), the set must be updated when new synthesis pathways are discovered. In contrast, we have preferred to use statistics extracted from training drug databases. This allows us to generate realistic sets of ligands whose properties will match the ones from the training database and which can be more easily synthesized. The program that is described in this manuscript can create completely de novo ligands using a growing approach (even when no actual ligands are known since only the coordinates of the active site are needed). It is also possible to start from a seed when prior experimental data are known. Depending on the size of the seed, this workflow can mimic a R-group strategy where the aim is to decorate a first hit, or a fragment-based drug discovery strategy where a small binding fragment is evolved to a small molecule. The approach used by the program resembles the one in SMoG in that it grows molecules in the active site of a protein by connecting small organic fragments [START_REF] Dewitte | SMoG: de Novo Design Method Based on Simple, Fast, and Accurate Free Energy Estimates. 1. Methodology and Supporting Evidence[END_REF] . However, SMoG ligands often lacked "chemical realism" and we aim to bring this "chemical realism" to our growth method. Since SMoG was released almost two decades ago, we decided to start the development of the new tool from scratch.

Moreover

Our program (OpenGrowth) is free, open-source, flexible and modular, easy to download, easy to install and easy to use. Amongst the important specifications we wanted, one can cite: taking into account protein flexibility from the beginning; having a flexible list of fragments; producing new ligands that are easy to synthesize and have good pharmacokinetic properties (ADME). To achieve the latter and bring the missing "chemical realism", it is crucial to choose new fragments in a reasonable way. We rely in the use of statistics of connection between fragments as defined in FOG [START_REF] Kutchukian | FOG: Fragment Optimized Growth algorithm for the de novo generation of molecules occupying druglike chemical space[END_REF] ; if the training drug database is big enough, significant information can be extracted and will produce compounds occupying a desired chemical space (i.e. drug-like molecules). Thus, our program is the first one to grow three-dimensional ligands using statistical analysis of training databases and having the possibility to take into account protein flexibility or to discover multi-target drugs, while being free, open-source and easy to use. The possibility to support fragment-based drug discovery with a seed approach makes OpenGrowth a very versatile tool for computational de novo drug design and open wide opportunities regarding its use.

Methods

Growing molecules

The growth begins by placing one of the fragments in the active site of the target protein; new fragments are then added to the current ligand until a given threshold is reached ("fragment" is used as a synonymous as "building block" throughout the text). Potential thresholds include molecular weight, number of atoms, or number of fragments (see Figure 1). The ligand is saved and then a new growth starts iteratively. The aim of the program is to grow hundreds of thousands of new ligands that can then be sorted by applying additional search criteria. The connection between the current state of the ligand and the next fragment is made by removing one hydrogen atom from each component and creating a new single bond between the two atoms to which the hydrogens were connected. The distance at which the new moiety is added depends on the type of the two atoms (for example, 1.54 Å between two sp [START_REF] Norinder | Prediction of ADMET Properties[END_REF] carbon atoms and 1.46 Å between sp 3 and sp carbon atoms). For each new fragment added, a rotameric search is performed about the new bond and the orientation with the lowest binding score and without steric clash is kept. A geometry optimization in the (fixed) active site of the protein can also be performed to relax the ligand structure. More details on the beginning of the growth, adding new fragments and optimizing the ligand geometry and its position are given in the Supporting Information.

The new ligand score is then compared to the previous score (before fragment addition), and a Monte-Carlo approach is used to accept or reject the fragment [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF] . If the new score is lower than the previous one, the new fragment is accepted; if the new score is higher than the previous one, the new fragment is accepted with a Metropolis criterion probability exp(-S/RT) where S=Snew-Sold is the difference of score and RT is a parameter given by the user. Thus, the growth tends to decrease the binding energy of the small molecules to the protein. Occasionally an increase in binding energy is allowed, which is useful in cases like crossing a tight channel. Figure 1. Principle of the OpenGrowth algorithm. After docking a primary fragment, new fragments are selected. For each one, once the best orientation is found a Metropolis criterion is used to either accept or reject the fragment. When the ligand reaches a given threshold property, it is saved, it can be filtered, and a new growth starts.

Drug library

Several aspects of the program use a drug library for training (for example, the construction of the fragment libraries or the evaluation of statistics of connection between fragments, see below).

Before going further, we want to describe quickly the library that was used. To prepare this set of molecules, we downloaded the full library of ChEMBL Drugstore (https://www.ebi.ac.uk/chembl/drugstore) current as of February 18 th 2014 (with the reference chembl_drugs-14_19-23-20) [START_REF] Bento | The ChEMBL bioactivity database: an update[END_REF] . This set contains 10,342 molecules, among which we kept the 8,145 containing a SMILES entry. This subset library was then cleaned by removing redundant molecules, single atom cations (Na + , K + , Mg 2+ , Ca 2+ , Li + , Zn 2+ , Al 3+ ), molecules with Si, Al, As, Pt, Bi, Sb atoms, molecules labelled as "inorganic" or "polymer", and small molecules such as bromoform, chloroform, carbon tetrachloride, etc… At the end, we obtained a library of 8,018 molecules. Not all these molecules have been approved by regulation agencies and some are still in clinical trials. However, we do not intend to capture by our algorithms the reasons why molecules may fail in clinical phase: our goal is first to create good protein binders. Since our approach uses statistical analysis of the drug library, we envisioned that statistics that are more accurate would be obtained with a larger library and we kept molecules that are still in clinical trials.

Fragment library

The choice of a proper fragment library is crucial for the success of a growth algorithm. As we will present later, we have developed a graphical interface for the user to choose which fragments he/she wants to use amongst a predefined list. To prepare this list, we sorted fragments in different classes.

First, we present the list of non-ring fragments in Figure 2. As a general principle, we tried to avoid denaturalizing the fragments. For example, an amide could theoretically be made during the growth by joining aldehyde and amine. However, we considered the nitrogen of an amide too different from that of an amine; thus, we added an amide fragment to the non-ring set. This approach is similar to the one used in FOG [START_REF] Kutchukian | FOG: Fragment Optimized Growth algorithm for the de novo generation of molecules occupying druglike chemical space[END_REF] . Some fragments presented in Figure 2 are trimers of other fragments (last line, for example urea could be made by joining amine, aldehyde and amine); the reason for the existence of the trimer fragments will be explained later. Some of the fragments in Figure 2 do not seem to be drug-like, such as carbonic acid (last line, first fragment); nonetheless, this fragment was actually found in the drug database. However it was found very rarely, thus it is very likely that the algorithm will almost never select it. A user can easily decide to remove this fragment from the list as explained in the user manual. Some organic fragments are not in the list since they are rarely found in drugs: azide (H-N3) was found 7 times in our drug library, sulfoxyde (O=SH2) was found 5 times and thionitrosyl (HN=S) only 2 times.

Moreover, we did not keep any fragments involving a phosphorus atom since most of these are unstable to hydrolysis; thus, a molecule containing a fragment with a phosphorus atom will probably not be considered drug-like at first (even if some drugs contain phosphorus atoms). One can manually enter new fragments in the list as well; we describe this process in the user manual. Nevertheless, the included set of non-ring fragments covers almost all the possibilities offered by organic chemistry. The choice of ring fragments is more complicated since many possibilities exist (note that in our algorithms two non-ring fragments cannot join to form a ring). We first extracted all the rings from the drug library by cutting all single bonds between rings and side chains. We achieved that by using the fragmenter tool from JChem 6.1.4 (ChemAxon, 2013, http://www.chemaxon.com). We can then assume that the rings that were found in our drug library (made of 8,018 molecules from ChEMBL) represent all the rings that can be found in drugs. We then classified this set of rings according to the number of rings fused together per fragment. We found the following numbers of different rings: mono-ring: 280; 2fused-rings: 444; 3-fused-rings: 361; 4-fused-rings: 171; 5-fused-rings: 77; 6-fused-rings: 28; 7-fusedrings: 8; 8-fused-rings: 6; 9-fused-rings: 5; 11-fused-rings: 3. We point out that in our definition, all double bonds connected to a ring are a part of the ring fragment (since we have cut single bonds to prepare the set of rings and we create single bonds during the growth). Thus, cyclopentane and cyclopentanone are considered different rings and are treated as such. To achieve this difference of treatment, we describe rings as SMARTS patterns instead of SMILES strings; SMARTS patterns allow us to add more information in the description of each atom. We created all the SMARTS patterns for the mono-rings and 2-fused-rings and validated them (for example, checking that a search for cyclopentane yielded no results that were actually cyclopentanone rings). Larger fragments (with three or more fusedrings) are not currently included. During the development of the program, we used only 21 mono-rings and six of the 2-fused-rings (in addition to the 24 non-rings fragments). These fragments were selected by prevalence in the drug library, and we kept those found the most often. This corresponds to a total of 113 fragments in the working library: indeed, since we use connection statistics between fragments (vide infra), a fragment such as pyridine must be described three times in the final library (one for each type of hydrogen). The 21 mono-rings and six 2-fused-rings which were used are displayed in Figure 3 and Figure 4. Expectedly, the most frequently occurring mono-ring and 2-fused-rings fragments are the same as the one found by Aldehi et al. and Taylor et al. [START_REF] Aldeghi | Two-and Three-dimensional Rings in Drugs[END_REF] (although not in the same order). Amongst the 2-fusedrings, one can find purine, xanthine and a -lactam ring constituting the core of cephalosporin antibiotics (closely related to penicillin). Should a user decide not to use the -lactam ring (e.g.) because the target protein is not related to antibiotics (whereas this fragment is mainly used in antibiotics), it can easily be discarded by not selecting it with the GUI. In addition to the fragments presented in Figure 3 and Figure 4, we have prepared all mono-rings and 2-fused rings fragment found more than 10 times in the drug library. This corresponds to a total of 58 mono-rings and 41 2-fused rings (413 different fragments in total). The full list of currently available fragments is displayed in the user manual. Moreover, the fragment preparation method is explained in the manual, allowing users to add new fragments according to their needs (3-fused-rings for example). The three-dimensional structures of each fragment were prepared by optimizing them with the Gaussian09 software (Rev. A.02) 37 at the M06-2X/aug-cc-pVTZ level of theory [START_REF] Zhao | The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals[END_REF] with an implicit description of water (IEF-PCM as implemented in Gaussian09) [START_REF] Tomasi | The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level[END_REF] . When it was relevant, we started from different conformations and kept the one with the lowest energy.

Connecting fragments

One of the lessons learned from SMoG was that fragments should not be randomly connected.

When no rules are used for selecting fragments, the output molecules can be unstable, toxic, nonsoluble and are also very likely difficult to synthesize. In 2009, our lab proposed to use chemoinformatic tools to analyze databases and to create de novo molecules occupying a given chemical space [START_REF] Kutchukian | FOG: Fragment Optimized Growth algorithm for the de novo generation of molecules occupying druglike chemical space[END_REF] , resulting in the so-called FOG algorithm. This approach has been used and improved here to fit the OpenGrowth needs. Briefly, we counted the connection between each pair of fragments in the drug database from ChEMBL. The probability that fragment i is connected to fragment j is then defined as:

𝑃 𝑖-𝑗 = 𝑁 𝑖-𝑗 ∑ 𝑁 𝑖-𝑓 𝑓
where Ni-j is the number of connections between fragments i and j and the sum in the denominator is made on all the fragments (more details are given in Supporting Information). In the original FOG algorithm, corrections had to be applied to the counting (see SI of [START_REF] Kutchukian | FOG: Fragment Optimized Growth algorithm for the de novo generation of molecules occupying druglike chemical space[END_REF]). However, the use of complete and more precise SMARTS patterns allows us to overcome this constraint. For example, amides are now

defined as [C;!R;X3;!$(C(=O)(N)N)](=[O;X1])-[N;!R;X3;!$(N(C(=O))C(=O));!$(N=C);!$(N(S(=O)(=O))C(=O))] (instead of [C&!R&X3](=O)-[N&!R&X3] in FOG).
In this new pattern, it is explicitly mentioned that the carbon is not part of a urea moiety (C(=O)(N)N), and that the nitrogen is not part of a N(C(=O))C(=O) moiety, an imine (N=C), nor a N(S(=O)(=O))C(=O) moiety. Since this amide fragment cannot be part of a urea group, it is necessary to define urea as a new fragment, explaining why some fragments in Figure 2 (last line) are trimers of other fragments. We can then create an N*N matrix with all the connections probabilities (where N is the number of strictly different fragments, N=113 for the selected fragments in this manuscript). Analysis of this matrix teaches us, for example, that a phenyl has a 29.3% chance of being connected to a methyl group, whereas a methyl has only a 6.4% chance of being connected to a phenyl ring. During the growth, new fragments are chosen according to these probabilities.

Consequently, the produced molecules will statistically "look like" molecules from the initial drug database. We will show later that the synthetic accessibility (as predicted by Sylvia) and the "drug-like"

properties (as predicted by QikProp) are subsequently superior to the ones from molecules randomly grown.

Searching for 3mers

To further increase the synthetic accessibility of new molecules, we implemented a 3mer-screen algorithm that was developed at the same time as FOG [START_REF] Kutchukian | FOG: Fragment Optimized Growth algorithm for the de novo generation of molecules occupying druglike chemical space[END_REF] . As we show further, this 3mer-screen also increases the "drug-like" properties of the new ligands (as predicted by QikProp). The algorithm is implemented as a sorting stage that is applied to grown ligands. The FOG growth is a two-state Markov chain that uses no information about connections to other fragments. Therefore, it is possible to connect three fragments together that are unlikely to be observed in that connected form in drugs, either because of synthetic difficulties, stability, or other reasons. For example, a -CH2-CH2-Br moiety is never found in our drug library. The 3mer-screen consists of first creating all the possible 3mers from the fragment library; the drug library is then screened for those 3mers, and we retain in a list the ones that are never found. The user can then manually add 3mers, such as acetals or ketals, which are known to be unstable. The result is a list of "forbidden 3mers", and when a newly grown molecule contains one of these 3mers the program warns the user about it. The principle of this screen is presented in Figure 5. A significant improvement has been made to the original algorithm: initially, rings were treated in a very general way with no specifications (SMARTS pattern [*&R]). We envisioned that this could lead to saving molecules that should have been discarded, especially concerning the connections around aromatic fragments. For example, it may be possible to find two phenyl substituents A and B in ortho to each other in drugs, but never in meta. We have thus used our SMARTS patterns to more accurately define 3mers involving rings, and the two patterns defining the case where the substituents are in ortho or in meta will now be seen as different cases. As a direct consequence, the list of forbidden 3mers is much longer; however, analyzing one new molecule against 420,000 forbidden 3mers in an AMD Opteron 6376 CPU (2.3GHz) takes only 20 seconds. Thus, searching for forbidden 3mers in the longer list does not significantly decrease the speed of the program.

Choosing fragments (GUI)

To facilitate the selection of fragments and the preparation of files needed by OpenGrowth (connection probabilities matrix, 3mer list, 3D structure files of the fragments…), we have developed a graphical user interface (GUI, see Figure 6). This interface is a program independent from OpenGrowth, and works under any version of Windows TM or Linux, using the Qt framework. It works in different stages: first, the user selects the "resource" and "output" directories and the database against which the probabilities will be calculated. We provide in the downloadable package our drug database (from ChEMBL) and the resource for the fragments (see user manual). Then, the user "updates" the database to count the occurrences of each fragment that are then displayed, as in Figure 6. Each fragment is a clickable button, allowing the user to select which one to keep; the "Fragment Files" button will create the list of the retained fragments. By default, non-ring fragments are always included and thus are not displayed on the GUI selection panel. If one wants to remove or add some of these fragments, he/she can manually alter the fragment list produced by the GUI. The "Probabilities" button reads the fragment list and creates the probability matrix, the "3-mer List" button creates the list of all the possible 3mer made by the fragment list and the "3-mers" button looks for the forbidden 3mers. All computational stages can be performed by the interface on a local machine or run on external servers. Figure 6. Graphical User Interface developed to select fragments, calculate connection probabilities, prepare the 3mer list and the 3D structure files.

Taking into account the protein flexibility

Using a single, rigid conformation of a protein is known to be a source of errors, leading to falsepositive and false-negative results [START_REF] Totrov | Flexible ligand docking to multiple receptor conformations: a practical alternative[END_REF] . To overcome this in the realm of docking, three approaches are mainly used 4b, 41 : (i) decreasing the van der Waals radii of atoms to give more space to the ligands ("softdocking") 42 ; (ii) sampling side chains of the protein residues, or removing them prior to adding them back after the ligand has been docked 43 ; (iii) docking the ligands in different conformations 40a, [START_REF] Zhao | FLIPDock: Docking flexible ligands into flexible receptors[END_REF] . We used an approach similar to that of the latter: in OpenGrowth, ligands can be grown simultaneously in several conformations of the protein. Since the program only needs .pdb files, the different conformations of the target protein can come from crystallographic structures of the same protein, structures from NMR experiments, or snapshots from molecular dynamics simulations. Moreover, all the input structures do not necessarily have to come from the same protein -the program can be used with different protein mutants, with the same protein from different bacteria or with strictly different targets (e.g.). At every stage of the growth, ligands orientations and geometries are optimized in each conformation independently (see Supporting Information). Therefore, to each protein conformation corresponds a specific ligand position, optimal for its own protein structure (see Figure 7). The score used for the Metropolis criterion is then either the lowest score in all structures, an arithmetic average or a

Boltzmann average (according to the user's choice):

𝑆 𝑎𝑣𝑒𝑟𝑎𝑔𝑒 = ∑ 𝑆 𝑖 * 𝑒 -𝛽𝑆 𝑖 𝑖 ∑ 𝑒 -𝛽𝑆 𝑖 𝑖 , 𝑆 𝑖 = 𝑆𝑐𝑜𝑟𝑒 𝑖𝑛 𝑐𝑜𝑛𝑓𝑜𝑟𝑚𝑎𝑡𝑖𝑜𝑛 𝑖
This latter approach shares similarities with previous work, such as those proposed by Fisher et al. where a Boltzmann average is used to take into account penalties from crystallographic occupancies [START_REF] Fischer | Incorporation of protein flexibility and conformational energy penalties in docking screens to improve ligand discovery[END_REF] or from Österberg et al. where Boltzmann averages are used to construct a grid [START_REF] Österberg | Automated docking to multiple target structures: Incorporation of protein mobility and structural water heterogeneity in AutoDock[END_REF] .

A modern view on protein flexibility consists of assuming that the protein "pre-exists in a number of energetically similar conformations" [START_REF] Teague | Implications of protein flexibility for drug discovery[END_REF] . When a ligand binds to a protein, it can bind to all these conformations. The most populated complex is not necessarily the one with the lowest binding free energy or the one from which the apo protein has the lowest folding free energy. The relevant energetic criterion is the sum of these two terms (plus the energetic deformation of the ligand, which is currently neglected). If we consider that the set of pre-existing conformations have a similar energy 47a , then the equilibrium between all the complexes relies only on the binding score, which justifies the use of a Boltzmann average. Hence, this is an initial way to take induced fit [START_REF] Koshland | Application of a Theory of Enzyme Specificity to Protein Synthesis[END_REF] into account, which supposes that the protein conformation will adapt around the ligand to lower the total Gibbs free energy. In our description, when the ligand is in a favorable position a higher weight is given to its score. Therefore, this specific protein conformation is considered to be adapted around the ligand as in the induced fit.

Scoring the protein/ligand binding energy

The current version of the program uses the same scoring function for both the rotameric search and the evaluation of the binding score for the Metropolis criterion. Two options are currently available:

the SMOG2001 scoring function [START_REF] Ishchenko | SMall Molecule Growth 2001 (SMoG2001): An Improved Knowledge-Based Scoring Function for Protein-Ligand Interactions[END_REF] and a more recent one called SMOG2015. More details on these functions are given in the Supporting Information. Since SMOG2015 is still in development, most of the results presented below were obtained with the SMOG2001 function (unless otherwise noted). Other scoring functions will be included in our program in future releases. However, since OpenGrowth is a modular platform and an open-source project, external contributors could easily implement their own scoring functions in the program (how to do so is explained in the user manual).

Computing accurate protein/ligand binding energy

Once OpenGrowth has generated in the order of 100,000 new ligands, it is necessary to screen these ligands with methods that are more accurate than scoring functions. Our test protein has been the HIV-1 protease, and we have used an approach proposed by Wright et al. to reproduce experimental binding Gibbs free energies [START_REF] Wright | Computing Clinically Relevant Binding Free Energies of HIV-1 Protease Inhibitors[END_REF] . We performed 50 independent molecular dynamics simulations of 5 nanoseconds for each protein-ligand complex, differing only by the initial velocities. Simulations were performed on GPU cards with Gromacs v5.0.2 [START_REF] Van Der | GROMACS: fast, flexible, and free[END_REF] , the Amber03 force field [START_REF] Duan | A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations[END_REF] and the TIP3P water model [START_REF]Transferable Intermolecular Potential Functions for Water, Alcohols, and Ethers. Application to Liquid Water[END_REF] .

Proteins were pre-processed by the Protein Preparation Wizard tool of Maestro (to add hydrogens and missing atoms and to minimize the energy) [START_REF]Protein Preparation Wizard, Schrödinger Release 2014-3[END_REF] . The correct protonation state of ligands was determined using Chimera 1.9 [START_REF] Pettersen | UCSF Chimera--a visualization system for exploratory research and analysis[END_REF] ; they were then optimized at the M06-2X/6-31+G** level of theory with an implicit description of water (IEF-PCM) using Gaussian09. Finally, RESP charges were computed at the HF/6-31G* level of theory with Gaussian09, and acpype (an interface to antechamber [START_REF] Wang | Automatic atom type and bond type perception in molecular mechanical calculations[END_REF] ) was used to generate the ligand topology using the General Amber Force Field 57 . The system was then solvated in a rhombic dodecahedron box with a distance of at least 11 Å between the complex and the box edges. A 10 Å cutoff was used for van der Waals interactions, whereas electrostatic ones were treated with PME [START_REF] Darden | Particle mesh Ewald: An N.log(N) method for Ewald sums in large systems[END_REF] with a similar "cut-off" of 10 Å. The energy of the solvated system was minimized until the maximum force reached 5.0 kJ.mol -1 .nm -1 , first with steepest descent and then with conjugate gradient. A primary NVT simulation of 200ps (dt=1fs, velocity-rescale thermostat [START_REF] Bussi | Isothermal-isobaric molecular dynamics using stochastic velocity rescaling[END_REF] , T=100K) was performed with the protein atoms being fixed (allowing the solvent to equilibrate). Then, another NVT simulation (without restraints)

was done with simulated annealing, the temperature going from 100K to 300K linearly in 400ps (dt=1fs, velocity-rescale thermostat); the system then evolved freely 100ps at 300K. A 500ps NPT simulation was subsequently carried out to equilibrate the pressure (dt=1fs, Nosé-Hoover thermostat 60 , Berendsen barostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] ). Finally, the production simulation was done (dt=2fs, Nosé-Hoover thermostat, Parinello-Rahman barostat [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] ). Bond lengths were constrained during the simulations using the LINCS algorithm [START_REF] Hess | P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation[END_REF] .

In the last nanosecond of each production simulation, 100 snapshots were extracted (every 10ps). The Amber12 package [START_REF] Case | The Amber biomolecular simulation programs[END_REF] was then used to compute the Molecular Mechanics/Poisson-Boltzmann Surface Area binding energy (MMPBSA) for the 5000 snapshots and the entropy contribution for 1000 snapshots (one every 50ps), using a single-trajectory approximation [START_REF] Hou | Molecular Dynamics and Free Energy Studies on the Wild-type and Double Mutant HIV-1 Protease Complexed with Amprenavir and Two Amprenavir-Related Inhibitors: Mechanism for Binding and Drug Resistance[END_REF] ; the results are then averaged. The ionic strength was fixed at 0.0M, whereas the surface tension and offset of the cavity of the Poisson-Boltzmann calculations were set to 0.0072 kcal/mol/Å 2 and 0.00 kcal/mol respectively [START_REF] Hou | Molecular Dynamics and Free Energy Studies on the Wild-type and Double Mutant HIV-1 Protease Complexed with Amprenavir and Two Amprenavir-Related Inhibitors: Mechanism for Binding and Drug Resistance[END_REF][START_REF] Hou | Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations[END_REF] (referred as Set1 below). We also performed calculations with surface tension and offset of the cavity set to 0.00542 kcal/mol/Å 2 and 0.92 kcal/mol respectively for comparisons (referred as Set2 below). The nonpolar optimization method of Amber12 ("inp" keyword) was set to 1; the use of the default value of 2 did not allow us to reproduce experimental results. For each complex, we verified that the distribution of the MMPBSA energies followed a normal distribution, ensuring that the 5000 snapshots represented a system at thermodynamic equilibrium. Each simulation took on average 6 hours to be completed (equilibration + production) whereas MMPBSA and entropy contributions took on average 8 hours on modern clusters. Note that since Gromacs and AmberTools are open-source softwares, all the simulations and analysis can be performed free of charge. The scripts and input files used for these simulations are available for download on the authors' websites and on the OpenGrowth webpage.

To validate this method, we applied it to a set of complexes with known crystallographic structures. Similarly to the work of Wright et al., we selected the 9 FDA-approved HIV-1 protease inhibitors [START_REF] Wright | Computing Clinically Relevant Binding Free Energies of HIV-1 Protease Inhibitors[END_REF] . To choose the complexes, we targeted sequences of the protein as close as possible to the consensus wild-type HXB2 [START_REF] Ratner | Complete nucleotide sequence of the AIDS virus, HTLV-III[END_REF] (same sequence as the PDB ID 1HXB). The list of retained structures is shown in Table 1. For some complexes, two orientations of the ligand (with occupancy smaller than 1.0) are given in the pdb structure. In such cases, simulations were started from both orientations and we kept the one with the lowest energetic values. The active site of the HIV-1 protease is made of two aspartate residues (ASP25) and there is a debate surrounding the protonation state of these amino-acids since four states are possible: the two residues from each chain can be both non-protonated, both protonated, or one or the other can be protonated. Recent studies have shown that the latter possibility is the most likely (see [START_REF] Mcgee | pH-REMD Simulations Indicate That the Catalytic Aspartates of HIV-1 Protease Exist Primarily in a Monoprotonated State[END_REF] and references cited therein). Thus, for each complex we have simulated the two possibilities (ASP25 of chain A protonated and not ASP25 of chain B, or vice versa). This set represents simulations of 15 complexes with two protonation states, corresponding to 1500 simulations of 5ns. The correlation coefficients (R) between the computational and experimental binding Gibbs free energies for different methods are presented in Table 2. We observed that for tight binding ligands (with Kd in the low picomolar (pM) range), the MMPBSA approach could fail to provide realistic estimations of the binding Gibbs free energy. This is the case when lopinavir is simulated with the complex from the PDB ID 1MUI (data not shown here) for which Kd=0.0013nM [START_REF] Stoll | X-ray crystallographic structure of ABT-378 (Lopinavir) bound to HIV-1 protease[END_REF] ; consequently, lopinavir was simulated with the complex from the PDB ID 2QHC instead (Kd=1.1nM) [START_REF] Šašková | Enzymatic and structural analysis of the I47A mutation contributing to the reduced susceptibility to HIV protease inhibitor lopinavir[END_REF] . Similarly, we observed that by removing darunavir (DRV) from the set of ligands used to benchmark our method (Kd=0.016nM with PDB ID 4HLA), the correlation was significantly increased. When Set2 of parameters are used for the MMPBSA calculations, the correlations between computations and experiments are slightly lower than with Set1 (see Table 2). The best correlation is obtained with Set1 of parameters, when entropy estimations are used, and when the lowest binging Gibbs free energy between the two protonation states is used (last column of the third line of Table 2, R=0.852, R 2 =0.726). For a majority of ligands (five of them), the lowest energy was coming from a structure where ASP25 from chain B is protonated (noted (3) in Table 2). Thus, it was the protonation state used for the protein when implemented in OpenGrowth in the future sections. The correlation between the computational (with Set1) and experimental binding Gibbs free energies for the eight kept FDA-approved HIV-1 protease inhibitors is presented in Figure 8. In the following, the "predicted dissociation constant" refers to the dissociation constant (Kd=exp(Gbinding/RT)) obtained using the correlation between computation and experiments, as shown in Figure 8.

Ligand

Technical details

OpenGrowth is written in C++ (6500 lines of code), and uses the OpenBabel library extensively [START_REF] O'boyle | Open Babel: An open chemical toolbox[END_REF] . It can typically grow 100 ligands/day/core in an AMD Opteron 6376 CPU (2.3GHz) and 170 ligands/day/core in an Intel Xeon E5-2640 v2 CPU (2.0GHz). The most computationally demanding part of the program is the geometry optimization of the ligands (see Supporting Information), and by tweaking the parameters (for example, performing less minimization steps) the speed can be increased.

OpenGrowth has been prepared in a modular manner, and each function can be updated to fit the user needs. We release it as an open-source project under the GNU GPL license. It is freely downloadable from the author websites and from http://opengrowth.sourceforge.net. OpenGrowth is easy to use ("./OpenGrowth.exe Input.inp") and uses a keyword-oriented input file. None of the parameters are hard-coded to allow the user to have more flexibility, and default values are proposed to simplify the process for new users. More details and explanations are provided in the Supporting Information of the current manuscript and on the program manual downloadable from the websites.

Results

Ability of the growth algorithm to reproduce transition probabilities

We wanted to verify that our program can grow new molecules with properties as desired, that is: connection probabilities between fragments, synthetic accessibility and pharmacokinetic properties, (d) Random connections (with OpenGrowth). Figure 9. Probabilities of connections from de novo molecules (grown with FOG2.0 or with OpenGrowth) versus the same probabilities in the drug database. Each point represents one probability, for example the probability that a phenyl ring is connected to an amine. We have underlined in purple the outliers discussed in the text.

Properties of output molecules

To evaluate the properties of the new molecules, different sets of ligands were grown in one conformation of the HIV-1 protease (coming from PDB ID 2AQU). The first set made with random connections between fragments; the second set made of molecules grown with FOG probabilities which have failed the 3mer screen; the third set also made with FOG probabilities but with molecules that passed the 3mer screen. During the growth, the program may generate several times the same molecule: 10% of the molecules are redundant with some others when random connections are used, whereas only 3% are redundant with others with FOG probabilities. After removing these redundant molecules, the three sets have sizes of 117,647, 53,759 and 21,754 molecules respectively (29% of all the molecules grown with FOG probabilities passed the 3mer screen). We compared the properties of these three sets with two drug libraries: the one that we have presented before (named ChEMBL drugs) and one made by Hutter with 2665 non-redundant molecules which "were chosen in a way that reflects the proportional distribution of all agents among the therapeutic categories" [START_REF] Hutter | Separating Drugs from Nondrugs: A Statistical Approach Using Atom Pair Distributions[END_REF] .

We first investigated the synthetic accessibility of the molecules by using the Sylvia program from Molecular Networks GmbH [START_REF] Boda | Structure and reaction based evaluation of synthetic accessibility[END_REF] . Sylvia returns a number between 1 and 10 after analyzing each molecule in terms of properties such as number of rings, number of stereocenters, similarities to a starting material database or molecular graph. The lower the score, the easier the synthesis of the molecule should be. We report Figure 10-(a) the distribution of Sylvia scores for the five sets of molecules. Molecules that were grown randomly (in yellow pentagones) have a high Sylvia score with a peak around 7.5. When FOG probabilities are used but the 3mer test failed (in orange circles), the distribution is shifted toward lower values of the score (peak at 6.5); passing the 3mer screen further decreases the score (in red triangles, peak at 5.5). For the two set of authentic molecules (in blue diamond and green squares), the peak is centered at 4.0. Thus, we validate here our strategy (connection probabilities and 3mer screen) to increase the synthetic accessibility of de novo ligands.

While this result was previously demonstrated for molecules grown in 2D with FOG [START_REF] Kutchukian | FOG: Fragment Optimized Growth algorithm for the de novo generation of molecules occupying druglike chemical space[END_REF] , we show here that we can access compounds with similar properties grown in the steric and energetic constraints of a binding pocket.

We then focused on the pharmacokinetic properties. Predicting ADME (Absorption, Distribution, Metabolization, Excretion) is a formidable challenge for computational drug design [START_REF] Norinder | Prediction of ADMET Properties[END_REF][START_REF] Wess | Medicinal Chemistry: Challenges and Opportunities[END_REF] ; the approach used by QikProp (Schrödinder, LLC) is to compute 43 descriptors or properties such as number of hydrogen bond donor atoms, polarizability, octanol/water partition coefficient or electron affinity [START_REF]QikProp, Small-Molecule Drug Discovery Suite 2014-3[END_REF] for each molecule. One must note that some key predictors of pharmacokinetic properties, such as metabolic liability and intrinsic clearance values, are not calculated by QikProp. Thus, even if QikProp cannot be used to "predict pharmacokinetic properties" (strictly speaking), it is a useful tool to identify "drug-like" molecules. For each property, 95% of drugs (from a pre-defined list) fall in a given range. Any time a computed property for a new molecule falls out of this range, one "star" is added. The number of computed stars for each molecule is then used as an estimation of its "drug-like" properties (the lower, the better). We present in Figure 10-(b) the analysis for the five sets of molecules. New ligands grown with FOG probabilities and which passed the 3mer screen (in red) are more likely to have a low number of stars than random molecules or molecules that did not pass the 3mer screen (yellow and orange). As expected, the two sets of drugs also have a higher probability of being given a low number of stars (blue and green). For instance, 66% of the ligands grown with FOG and which passed the 3mer screen have 0 or 1 stars, whereas this value falls to 53% for FOG molecules that failed the 3mer screen and 44% for molecules that were randomly grown. For the Hutter and ChEMBL databases, there are 83% and 79% molecules that have 0 or 1 stars. For the randomly grown library, 10% of the molecules have 6 stars or more, which falls to 7% for FOG molecules that failed the 3mer screen, 3% for FOG molecules that passed the 3mer screen and less than 2% for the two drug databases. So here again, the selection of molecules in the context of a binding pocket does not preclude the generation of drug-like compounds, and the implementation of connection probabilities and 3mer screen increases the chance to make new ligands with good ADME properties. It must be pointed out that even if the size of the five sets are not the same, we have checked that the same distribution (for Sylvia and for QikProp) is obtained when we randomly down-sample the libraries. We then analyzed the diversity of the grown ligands. Of the 21,754 molecules created with FOG probabilities that passed the 3mer screen, 3% were redundant with others in the library. We then computed the Tanimoto score between all pairs of the remaining 21,138 molecules (corresponding to more than 223 million of scores). Tanimoto score is a measure of similarity between two molecules: it is a number between 0 (very different) and 1 (very similar). The distribution of the observed scores is presented in Figure 11. The peak of the distribution is around 0.2, which signifies highly different molecules. Thus, even if the growth is biased towards a given chemical space, the OpenGrowth algorithm grows molecules that occupy a large part of this chemical space.

Predicting binding energy of de novo ligands

To evaluate if our program can grow molecules that are predicted to be good binders, we took the 75,000 ligands grown with FOG connection probabilities that were used in the previous section.

They were grown in one conformation of the HIV-1 protease (from PDB ID 2AQU) using a threshold of 20 fragments and the SMOG2001 scoring function. Since we wanted to focus here on the binding free energy, we also considered molecules that failed the 3mer-screen. Scores ranged from -90 to 0 in SMOG2001 internal units which correlate to binding free energy (for comparison, scores for the eight authentic HIV-1 protease inhibitors ranged from -44 to -24 SMOG2001 units). We picked 10 of the bests de novo ligands and applied the MD protocol defined previously (in section "Computing accurate protein/ligand binding energy"). We display in Table 3 some of the predicted properties for these molecules and the highest Tanimoto score between each ligand and the full ZINC library (which is made of 19 million purchasable compounds as updated on November 28 th 2014) [START_REF] Irwin | ZINC--a free database of commercially available compounds for virtual screening[END_REF] . The structures of the nearest neighbors are presented in Supporting Information: some similarities were observed after visual inspections, however no significant similarities could be found. The most similar molecule came from the analysis of Ligand014_18 for which ZINC67956366 was the closest. The four molecules with submicromolar predicted Kd are displayed in Figure 12. 3. Ten of the best de novo ligands (according to SMoG2001) grown. The SMOG2001 score is in internal units. MW is the molecular weight in g/mol. The 3mer column displays the success or failure of the 3mer screen. "Predicted Kd" is the predicted dissociation constant calculated using Figure 8. The Sylvia score is in internal units and ranges between 1 and 10. QikProp displays the number of "stars" per molecule. The Tanimoto column shows the highest score for the comparison between each molecule and the full ZINC library. Figure 12. Four predicted sub-micromolar ligands to HIV-1 protease.

Ligand

Therefore, in less than a week we were able to generate new molecules that are predicted to bind to the HIV-1 protease in the nanomolar range. These molecules might constitute interesting starting points for finding a first hit, especially Ligand014_279 and Ligand014_18 due to their low Sylvia scores for synthetic accessibility. We noted, however, that some of the top scored ligands (according to SMOG2001) are in fact predicted to be poor binders by MD simulations, with predicted Kd above 100M. This should certainly be attributed to the limitations of predicting power of SMOG2001, as will be discussed later.

Seeding a fragment

Most of drug discovery projects start with the identification of a hit that is later optimized to a lead. This is especially true in fragment-based drug discovery (which has emerged as a promising alternative to high-throughput screening) where low-molecular weight fragments that binds the target protein are identified and grown to larger molecules. A growth with OpenGrowth can start with a seed to mimic this kind of strategy. Thus, OpenGrowth can provide a support to fragment-based drug discovery projects by generating new ideas on how to evolve a fragment into a small molecule. During the development of atazanavir (sold under the name Reyataz by Bristol-Myers Squibb), a first scaffold was identified (see Figure 13) [START_REF] Fässler | Novel pseudosymmetric inhibitors of HIV-1 protease[END_REF] . Subsequent variations of this aza-dipeptide moiety and lead optimization led to the approval of the drug by the FDA in 2003 [START_REF] Fässler | Aza-Peptide Analogs as Potent Human Immunodeficiency Virus Type-1 Protease Inhibitors with Oral Bioavailability[END_REF] . We grew 170,000 ligands using the seed displayed in Figure 13. Growth was only allowed from the marked atoms that were identified in the input file as hydrogen atoms from an amide or a methyl, respectively. Thus, we utilized FOG probabilities to grow molecules. Ten of the best ligands are listed in Table 4 with some of the predicted properties and the highest Tanimoto scores with respect to the full ZINC library. Eight of them have sub-micromolar predicted Kd and could have been served as starting point for the optimization to provide new hits. Four of them are displayed in Figure 14. We note that starting from a seed allow obtaining more homogeneous predicted Kd: eight out of the 10 bests are in the hundreds of nanomolar range, and the two remainings are in the low micromolar range. When no seed was used (Table 3), predicted Kd are very diverse and range from 18nM to 110M. 4. Ten of the best de novo ligands (according to SMOG2001) grown around a seed. The SMOG2001 score is in internal units. MW is the molecular weight in g/mol. The 3mer column displays the success or failure of the 3mer screen. "Predicted Kd" is the predicted dissociation constant calculated using Figure 8.

Ligand

The Sylvia score is in internal units and ranges between 1 and 10. QikProp displays the number of "stars" per molecule. The Tanimoto column shows the highest score for the comparison between each molecule and the full ZINC library. Figure 14. Four predicted sub-micromolar ligands to HIV-1 protease grown around a seed.

We also compared each one of the 170,000 grown ligands with the atazanavir structure by computing the Tanimoto scores. The most similar molecule (Seed101_146) has a Tanimoto score with atazanavir of 0.80, a SMOG2001 score of -62.18 units and is ranked (according to SMOG2001 score) 28,713 in the list of all molecules. Its structure is displayed on Figure 15; we also show the superimposition of the positions of atazanavir and Seed101_146 in the active site. When comparing the two molecules (see atazanavir in Figure 13), several common features appear: (1) the phenyl-pyridine moiety was grown by OpenGrowth at the same position than the original structure. This moiety (present at one of the extremity of the active site) develops hydrophobic interactions with Val82 and Pro81 from chain A and Ile50 from chain B. (2) On the other side of the molecule, an isopropyl moiety was grown in place of a tertbutyl in atazanavir, the two being very similar. (3) An amine was added in beta of the amide on the left whereas an amide is present in atazanavir, the two moieties sharing the nitrogen atom at the same place. Overall, Seed101_146 shares several interesting features with the FDA-approved drug. However, it must be pointed out that Seed101_146 is predicted as a weak binder, which suggests room for improving our scoring method. The use of more accurate scoring functions than SMOG2001 will remove this limitation and open many interesting applications for this program.

Using the protein flexibility option

We first used the protein flexibility approach with 10 conformations of the HIV-1 protease and the SMOG2001 scoring function. The 10 conformations were obtained by performing 10 independent molecular dynamics simulations of the apo protein as explained previously (in section "Computing accurate protein/ligand binding energy") and starting from PDB ID 2AQU. We clustered the last four nanoseconds of each simulation and kept the central structure of each cluster (which is the conformation with the smallest average root-mean-square deviation to all the other structures in a given trajectory). We grew 10,000 molecules with OpenGrowth and picked 10 of the best structures to compute the binding free energy more accurately (see Table 5). However, we found that the ligands grown with this method were smaller than the ones previously presented: the lowest SMOG2001 score is only of -50.20 units and the molecular weight is on average 220 g/mol (whereas Table 3 and Table 4 show ligands of 540 g/mol). Consequently, the predicted Kd are higher (Kd=638nM for Flex073_23 or Kd=1.81M for Flex071_6). The reason for this behavior is very likely the hard-wall potential used with the SMOG2001 scoring function that imposes many more constraints; we envisioned that a scoring function lacking this constraint would probably perform better for predicting good binders. We computed the Tanimoto score between each one of these ligands and the full ZINC library and the structures of the nearest neighbors are presented in Supporting Information: similarities can be found, but none of the grown ligands are identical to molecules contained in the ZINC library (even for pairs with a Tanimoto score of 1.00). 5. Ten of the best de novo ligands (according to SMOG2001) grown in 10 structures of the protein.

Ligand

The SMOG2001 score is in internal units. MW is the molecular weight in g/mol. The 3mer column displays the success or failure of the 3mer screen. "Predicted Kd" is the predicted dissociation constant calculated using Figure 8. The Sylvia score is in internal units and ranges between 1 and 10. QikProp displays the number of "stars" per molecule. The Tanimoto column shows the highest score for the comparison between each molecule and the full ZINC library.

As stated previously, we are currently developing a new scoring function called SMOG2015.

Besides using a more accurate training set, another kind of atom typing and three shells instead of two, this scoring function uses a repulsive term instead of a hard-wall potential for protein-ligand interaction.

We present in Table 6 results of preliminary tests while growing compound in one and in five target structures. For both cases, we present the SMOG2015 score and the molecular weight of the 10 best ligands. It appears clearly that while using SMOG2015 with five structures, it is possible to grow ligands of the same size as with only one structures (on average, the molecular weight of the 10 best ligands is 598 g/mol when growing in one structure and 513 g/mol while growing in five structures). Thus, SMOG2015 can overcome the limitations of SMOG2001 while using the protein flexibility option. 6. Preliminary results obtained when growing ligands with the SMOG2015 scoring function in one and in five target structures. The SMOG2015 score is in internal units. When 5 structures were used, we have used an arithmetic average of the scores. MW is the molecular weight in g/mol.

Ligand

Discussion and Conclusions

We have presented here the foundations of a new lead design tool named OpenGrowth for use in the growth of ligands in the active site of a protein. Due to the open-source nature of this project, external contributors have the opportunity to improve the program. The efficiency of our algorithm (more than 150,000 new ligands can be grown in one week using 128 cores) allowed us to quickly identify new structures with a predicted Kd with the HIV-1 protease in the nanomolar range. This makes OpenGrowth a very fast and efficient tool for drug discovery, allowing one to discover novel scaffolds for further optimization. Ligands grown by our program must be analyzed and refined before being tested experimentally, but several of the structures presented here need only small modifications to be synthesized for trial. For example, one could find that the carboxylic moiety in Ligand014_279 (see Figure 12) is unsuitable for HIV-1 protease, and thus decide to remove it. Similarly, Seed042_1142 and Seed004_360 (see Figure 14) could be optimized to be more stable in vivo. Importantly, the possibility of growing molecules from a seed in OpenGrowth can dramatically speed up the search of variations around a first hit shortening the time needed for the discovery of leads or candidates. As such, it is a perfect in silico tool to tackle the problem of fragment-based drug discovery. In addition, this strategy can also be used to build ligands to metalloproteins: describing correctly the metal-ligand interactions with a scoring function is challenging, thereby in the past our group has used a seeding approach to discover new ligands to human carbonic anhydrase [START_REF] Grzybowski | Combinatorial computational method gives new picomolar ligands for a known enzyme[END_REF] .

The growth of molecules in OpenGrowth is biased to produce compounds that share properties with a training drug database, specifically the probability that two substructures are connected. To further increase the synthetic accessibility of the new molecules, we have also implemented the previously described 3mer screen. We have shown that when used together, these two options allow generating much better synthetic accessibility and pharmacokinetic/"drug-like" properties for the de novo molecules than for molecules created by random growth. We plan to include new options in OpenGrowth to compute a synthetic accessibility score or pharmacokinetic/"drug-like" properties directly. These properties could be used during the growth process using multi-objective optimization: for example, addition of a new fragment might not only have to satisfy the Metropolis criterion, but also not increase the accessibility score of the molecule past a certain threshold. If the score increases too much after the addition of three fragments (e.g.) then this set of three fragments might be rejected.

Such extended analysis will likely increase the success of the growth process.

One of the current limitations of OpenGrowth is the scoring function. Even though SMOG2001 has the same accuracy as the well-known GlideScore-SP with a correlation coefficient of R=0.45 with experimental data (see Supporting Information), this may be enough to be used to guide growth with the Metropolis criterion or to efficiently sort structures. Moreover, the hard-wall potential may be a source of rejection for too many fragments. To overcome this, we are currently developing a more accurate scoring function in the same spirit as SMOG2001 (knowledge-based potential with a multi-shell approach), called SMOG2015. Its preliminary version is already implemented in OpenGrowth, and full details will be provided in a future publication. Briefly, a repulsive term and ligand entropy contribution (by the mean of the number of rotatable bonds [START_REF] Huey | A semiempirical free energy force field with charge-based desolvation[END_REF] ) have been included in the new scoring function to increase its predictive power. Moreover, since many more complexes are available in the Protein Data Bank than in 2001, statistics that are more accurate have been obtained (see Supporting Information).

Preliminary results using this scoring function are promising, since it allowed us to grow relatively large compounds in five target structures, whereas SMOG2001 failed to do the same. Due to the open-source nature of OpenGrowth, the implementation of existing scoring functions in our program by external contributors is also possible (and is explained in the user manual). In the future, the rotameric search, the calculation of the score for the Metropolis criterion and the final energy prediction may use different energy functions. For example, the final energy prediction could be similar to a MMPBSA calculation.

Regarding the protein flexibility, future releases of OpenGrowth could include weights for the different protein structures according to their free energies; moreover, in case we do not want the ligand to bind to off-target proteins, it is also possible to include a penalty for those proteins.

Besides the demonstrated applications of OpenGrowth presented here, we want to underline that the fields of applications of this program are very wide. For example, peptides can theoretically be grown using the current version of OpenGrowth, but we plan to simplify this process by preparing the needed fragments; this approach could serve to design new inhibitors to proteases by including nonnatural amino-acids in the fragment library. One can also think of growing compounds to strengthen protein-protein interactions. Future releases of the program will allow the possibility to grow covalent ligands using the seeding approach, which would be a major breakthrough in the realm of drug design.

Thus, the ease of use, the flexibility and the open nature of the OpenGrowth program makes it a strong advancement over previous tools for drug design.

Supporting Information Available

More details on the growth (preparation and addition of a fragment), the optimization of the ligands, score computation, the FOG probabilities and the nearest neighbors to new ligands are available in the Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org.

Corresponding Author Information

Email: shakhnovich@chemistry.harvard.edu. Phone number: +1-617-495-4130.

Table of Contents Graphic

  , while the current version of LigBuilder can only work with one target, the program presented here has the possibility to take into account protein flexibility or to create multitarget compounds by growing ligands in different target proteins at the same time. Finally, since our program is open-source external contributors can provide new ideas or implement other scoring functions whereas the sources of LigBuilder are not accessible.

Figure 2 .

 2 Figure 2. List of non-ring fragments.

Figure 3 .

 3 Figure 3. List of mono-ring fragments used. The number of instances of each ring in the drug library is

Figure 4 .

 4 Figure 4. List of 2-fused-ring fragments used. The number of instances of each ring in the drug library is

Figure 5 .

 5 Figure 5. Principle of the 3mer-screen algorithm.

Figure 7 .

 7 Figure 7. Growing in five conformations of a protein (the five conformations were previously aligned).

Figure 8 .

 8 Figure 8. Correlation between computation and experiments for the binding Gibbs free energy for eight

  Figure 9-(b) and Figure 9-(d) and the non-perfect correlation in Figure 9-(a) and Figure 9-(c).

  (a) FOG probabilities (with FOG2.0). (b) Random connections (with FOG2.0). (c) FOG probabilities (with OpenGrowth).

  (a) Synthetic accessibility. (b) ADME properties (using QikProp).

Figure 10 .

 10 Figure 10. Properties of the output molecules.

Figure 11 .

 11 Figure 11. Distribution of the Tanimoto score between all pairs of 21,138 molecules. A score of 1.0

Figure 13 .

 13 Figure 13. Atazanavir and the seed used for growth.

  (a) Seed101_146. (b) Superimposition with the position of atazanavir.

Figure 15 .

 15 Figure 15. Most similar grown structure to atazanavir and superimposition of structures in the active

Table 1 .

 1 List of retained HIV-1 protease-inhibitor complexes to validate the computational estimation of the binging Gibbs free energy.

		Abbreviation	Mutations	PDB ID	Kd (nM)	Reference
	Amprenavir	APV	V3I	1HPV	0.6	69
	Atazanavir	AZV	V3I	2AQU	0.48	70
	Darunavir	DRV	V3I	4HLA	0.016	71
	Indinavir	IDV	V3I	1HSG	2.45	72
	Lopinavir	LPV	V3I/I47A	2QHC	1.1	73
	Nelfinavir	NFV	V3I	1OHR	2.0	74
	Ritonavir	RTV	V3I/S37N	1HXW	0.015	75
	Saquinavir	SQV	Wild-Type	1HXB	0.12	76
	Tipranavir	TPV	V3I/Q7K	2O4P	0.019	77

  Set1 refers to values of 0.0072 kcal/mol/Å 2 and 0.00 kcal/mol for the surface tension and offset of the cavity respectively. Set2 refers to values of 0.00542 kcal/mol/Å 2 and 0.92 kcal/mol.

		MMPBSA (2)	MMPBSA + Entropy (2)	MMPBSA (3)	MMPBSA + Entropy (3)	MMPBSA (Best)	MMPBSA + Entropy (Best)
	With 9 ligands (Set1)	0.483	0.563	0.695	0.598	0.597	0.582
	Without DRV (Set1)	0.695	0.822	0.835	0.766	0.808	0.852
	Without DRV (Set2)	0.696	0.820	0.823	0.749	0.805	0.841
	Table 2. Correlation coefficient (R) obtained with different methods. (2) means that ASP25 of chain A is
	protonated and not ASP25 from chain B. (3) means the opposite: ASP25 of chain B is protonated and not
	ASP25 of chain A. (Best) means that the lowest binding free energy is kept, coming either from (2) or
	from (3).

(10 structures) SMOG2001 MW 3mer Predicted Kd (nM) Sylvia QikProp Tanimoto

	Flex073_23	-50.20	246	Fail	638	4.14	0	0.84
	Flex095_24	-50.15	262	Fail	5,530	4.90	0	0.73
	Flex097_66	-44.66	234 Pass	4,940	4.38	2	0.72
	Flex112_6	-44.30	206 Pass	587,000	4.32	0	0.84
	Flex120_64	-44.19	276 Pass 81,200,000 4.75	3	0.82
	Flex097_42	-43.62	186 Pass	943,000	3.09	1	1.00
	Flex063_14	-43.43	193	Fail	2,093,000	3.77	1	0.89
	Flex088_16	-42.66	203 Pass	12,300	4.40	1	0.79
	Flex071_6	-42.49	197 Pass	1,810	4.96	2	1.00
	Flex100_7	-42.39	202 Pass	4,290	4.01	6	1.00
	Table

Acknowledgments

This work was supported by Defense Advanced Research Projects Agency (DARPA) contract HR0011-11-C-0093. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. The authors thank Dr. Peter Kutchukian for his comments on the manuscript. N.C. thanks the Région Rhône-Alpes for an Explora'Pro fellowship. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.


Nicolas Chéron

Chenchen Yu

Abimbola O Kolawole

Eugene I Shakhnovich

Christiane E Wobus
email: cwobus@umich.edu

Repurposing of rutin for the inhibition of norovirus replication

come

Main text

The development cost of new drugs has exploded over the last decades and is now estimated to be at least $2 billion per drug [27]. At the same time, the number of new molecular entities approved by FDA has saturated to about 25 per year [START_REF] Kinch | An overview of FDA-approved new molecular entities: 1827-2013[END_REF]. The success rate of new molecular entities is estimated to be 10% from Phase II to market, and 50% from Phase III to market. These values increase to 25% and 64% for repositioned drugs [START_REF] Reuters | Knowledge-based drug repositioning to drive R&D productivity[END_REF]. Therefore, an increased effort is being made towards drug repurposing [32, [START_REF] Sardana | Drug repositioning for orphan diseases[END_REF][START_REF] Xu | Database identifies FDA-approved drugs with potential to be repurposed for treatment of orphan diseases[END_REF], with demonstrated success in part due to the fact that different diseases can have related biological origin [START_REF] Reuters | Knowledge-based drug repositioning to drive R&D productivity[END_REF]37,46]. Since toxicology and pharmacokinetic data of repurposed molecules are already available, the development process is accelerated and the risk of failures is lowered.

Noroviruses (NoV) cause almost 80% of non-bacterial gastroenteritis in humans and spread quickly in confined environments, leading to 71,000 hospitalizations and 800 deaths every year in the USA [START_REF] Blanton | Molecular and epidemiologic trends of caliciviruses associated with outbreaks of acute gastroenteritis in the United States, 2000-2004[END_REF][START_REF] Hall | Norovirus disease in the United States[END_REF][START_REF] Patel | Systematic literature review of role of noroviruses in sporadic gastroenteritis[END_REF][START_REF] Scallan | Foodborne illness acquired in the United States--major pathogens[END_REF]. The impact is even greater in resource poor nations, where an estimated 200,000 children under 5 years of age die every year from NoV infections [START_REF] Patel | Systematic literature review of role of noroviruses in sporadic gastroenteritis[END_REF]. However, no specific treatments to control or prevent norovirus illness are approved by the FDA to date. Norovirus vaccine efforts are in phase I clinical trials [START_REF] Treanor | A Novel Intramuscular Bivalent Norovirus Virus-Like Particle Vaccine Candidate-Reactogenicity, Safety, and Immunogenicity in a Phase 1 Trial in Healthy Adults[END_REF], however vaccines cannot be used in some at risk populations, such as immunocompromised patients [17]. Moreover, any effective vaccine design needs to overcome the large genetic diversity and rapid mutation rate of HuNoVs [START_REF] Kroneman | Proposal for a unified norovirus nomenclature and genotyping[END_REF]. Thus, the development of low-molecular weight drugs with anti-norovirus activity has emerged as another strategy to treat or prevent norovirus infections. To date, three of the nine proteins encoded by the (+)-sense, singlestranded RNA genome have been researched extensively for their potential as antiviral targets [START_REF] Arias | Progress towards the prevention and treatment of norovirus infections[END_REF]: the capsid protein VP1 [START_REF] Feng | Library screen for inhibitors targeting norovirus binding to histo-blood group antigen receptors[END_REF][START_REF] Rademacher | Targeting norovirus infection-multivalent entry inhibitor design based on NMR experiments[END_REF], the protease [START_REF] Kim | Broad-spectrum antivirals against 3C or 3C-like proteases of picornaviruses, noroviruses, and coronaviruses[END_REF][START_REF] Mandadapu | Potent inhibition of norovirus 3CL protease by peptidyl alpha-ketoamides and alpha-ketoheterocycles[END_REF][START_REF] Tiew | Design, synthesis, and evaluation of inhibitors of Norwalk virus 3C protease[END_REF], and the RNA-dependent RNA-polymerase (RdRp) [START_REF] Costantini | Antiviral activity of nucleoside analogues against norovirus[END_REF][START_REF] Rocha-Pereira | E)-2-styrylchromones as potential anti-norovirus agents[END_REF][START_REF] Zamyatkin | Structural insights into mechanisms of catalysis and inhibition in Norwalk virus polymerase[END_REF]. In the absence of a robust HuNoV cell culture system [START_REF] Karst | Advances in norovirus biology[END_REF], murine NoV (MNV) has proved to be an attractive and widely used surrogate virus for HuNoV, including for drug development efforts [START_REF] Rocha-Pereira | The Viral Polymerase Inhibitor 2'-C-Methylcytidine Inhibits Norwalk Virus Replication and Protects against Norovirus-Induced Diarrhea and Mortality in a Mouse Model[END_REF].

To identify molecules that potentially could be repurposed as a broad spectrum norovirus antiviral, we chose molecules that have already been used in humans and docked them in silico into the MNV RdRp protein (PDB ID: 3UPF [START_REF] Mastrangelo | Structure-Based Inhibition of Norovirus RNA-Dependent RNA Polymerases[END_REF]). The "ZINC In Man" library was downloaded in SMILES format (http://zinc.docking.org/browse/subsets/special) [START_REF] Irwin | ZINC--a free database of commercially available compounds for virtual screening[END_REF][START_REF] Irwin | ZINC: a free tool to discover chemistry for biology[END_REF]. The 11421 molecules were converted to 15990 3D structures with LigPrep [START_REF] Butler | Modulation of dendritic cell phenotype and function in an in vitro model of the intestinal epithelium[END_REF] by preparing different stereoisomers or protonation states. This in silico library was first docked with Glide-SP on chain B of the MNV RdRp [2, [START_REF] Friesner | Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy[END_REF][START_REF] Halgren | Glide: a new approach for rapid, accurate docking and scoring. 2. Enrichment factors in database screening[END_REF]. The 2500 best molecules (i.e. with the lowest scores) were kept and docked with the extraprecision mode (XP) of Glide [START_REF] Friesner | Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes[END_REF] on the same grid. The 37 molecules, which can be represented as 90 states (stereoisomers or tautomers), was next docked on the three chains of the MNV RdRp with the QM-Polarized Ligand Docking approach and the scores were averaged [START_REF]QM-Polarized Ligand Docking protocol[END_REF][START_REF] Cho | Importance of accurate charges in molecular docking: quantum mechanical/molecular mechanical (QM/MM) approach[END_REF]. This set of 90 molecules was also docked into the human norovirus RdRp (PDB ID: 1SH2 [START_REF] Ng | Crystal structure of norwalk virus polymerase reveals the carboxyl terminus in the active site cleft[END_REF]). The docking outcome was then analyzed, and we removed the molecules which appeared not promising for testing as a potential drug due to their ubiquity in cells (for example NADH or sugars). Sixteen molecules with the lowest score for MNV RdRp and the accompanying scores for HuNoV RdRp are presented in Table 1.

Based on the scores for murine and human NoV RdRp and the price of the molecules, four compounds were selected for experimental testing in cell culture: ZINC28639308 (AD132644), ZINC11726792 (Sinefungin), ZINC72206342 (Rutin), ZINC04216682 (Kasugamycin). The predicted dissociation constants of these molecules with the MNV RdRp are 1.4, 66, 78 and 588 nM, respectively. Rutin hydrate (Sigma, Cat. No. R5143) and kasugamycin hydrochloride from streptomyces kasugaensis (Sigma, Cat. No. K4013) were freshly prepared in DMSO as a 10 mM stock solution or in water as a 1 mM stock solution, respectively. Sinefungin (Sigma, Cat. No. S8559) and AD132644 (AnalytiCon Discovery, Cat. No. NP-005514) were dissolved in water as a 26.22 mM stock solution, or in DMSO as a 10 mM stock solution, respectively, and stored at 4 o C. Since no reproducible HuNoV cell culture model was available [START_REF] Duizer | Laboratory efforts to cultivate noroviruses[END_REF][START_REF] Karst | Advances in norovirus biology[END_REF], we tested the antiviral activity of these compounds against MNV in the murine macrophage cell line RAW 267.4 (ATCC; Manassas, VA). Cells were maintained as described previously [START_REF] Wobus | Replication of Norovirus in cell culture reveals a tropism for dendritic cells and macrophages[END_REF], and the plaque-purified MNV-1 clone (GV/MNV1/2002/USA) MNV-1.CW3 (referred herein as MNV-1) was used at passage 6. Our first goal was to determine the cellular toxicity of the selected compounds in RAW 264.7 cells using the WST-1 assay (Roche, San Francisco, CA), which measures mitochondrial dehydrogenase activity [16]. Concentrations that maintained cell viability at or above 80% were considered non-toxic. Concentrations equal or above 100 M for rutin (Figure 1a), 0.5 M for kasugamycin (Figure 1b), 1 mM for sinefungin (Figure 1c), and 5 M for AD132644 (Figure 1d) maintain cellular viability.

In order to identify whether the selected drugs exhibited antiviral efficacy against NoV, RAW 264.7 cells were infected with MNV-1 (MOI of 2.0) in the presence of rutin, kasugamycin, sinefungin, and AD132644 at previously determined non-toxic concentrations, or equal amounts of vehicle control (i.e., H2O or DMSO) (Figure 2). RAW 264.7 cells were plated at 5 x 10 5 per well in 12-well plates the day before. Infection occurred for 1 h on ice, after which the viral inoculum was removed by washing cells with ice-cold PBS three times before adding media containing one of the four drugs at indicated concentrations. Eight hours later, viral titers were determined by plaque assay as described previously [START_REF] Gonzalez-Hernandez | Chemical Derivatives of a Small Molecule Deubiquitinase Inhibitor Have Antiviral Activity against Several RNA Viruses[END_REF]. A statistically significant difference between rutin and vehicle-treated cells was observed at 1 M and 100 M of drug (Figure 2a). However, no antiviral effect was observed at any of the concentrations tested for kasugamycin (Figure 2b), sinefungin (Figure 2c), and AD132644 (Figure 2d). Taken together, these data suggest that rutin but not kasugamycin, sinefungin and AD132644 have antiviral efficacy against MNV-1 in cell culture.

Rutin is a flavonoid which belongs to the family of vitamin C2 and is found in many vegetables (olive, asparagus), fruits (raspberry, plum) or cereals (buckwheat) (http://phenolexplorer.eu/contents/polyphenol/296). It is a compound with known antioxidant [START_REF] Yanga | In vitro antioxidant properties of rutin[END_REF], antiinflammatory [START_REF] Guardia | Anti-inflammatory properties of plant flavonoids. Effects of rutin, quercetin and hesperidin on adjuvant arthritis in rat[END_REF] and anti-cancer effects [START_REF] Deschner | The effect of dietary quercetin and rutin on AOM-induced acute colonic epithelial abnormalities in mice fed a high-fat diet[END_REF][START_REF] Lin | Rutin inhibits the proliferation of murine leukemia WEHI-3 cells in vivo and promotes immune response in vivo[END_REF] and some of its derivatives are commercialized under the name Venoruton. To get more insight on the mechanism of action of rutin, we performed RNA synthesis assays using recombinantly expressed and purified MNV-1 and GII.4 HuNoV RdRps [START_REF] Yunus | The murine norovirus core subgenomic RNA promoter consists of a stable stem-loop that can direct accurate initiation of RNA synthesis[END_REF].

Surprisingly, rutin exhibited no inhibitory activity against either RdRp at all concentrations tested up to 5 M (see Supporting Information). In the same assay, benzodiathiadine, a potent non-nucleoside inhibitor of the HCV RdRp [START_REF] Dhanak | Identification and biological characterization of heterocyclic inhibitors of the hepatitis C virus RNA-dependent RNA polymerase[END_REF], inhibited RNA synthesis by the HCV RdRp at low micromolar concentrations (data not shown). To determine whether rutin acts as a prodrug, its activity was tested in a mammalian cell-based assay where transiently expressed RdRp stimulates retinoic acid-inducible gene I (RIG-I)-dependent reporter luciferase production via the beta interferon promoter [START_REF] Subba-Reddy | VPg-primed RNA synthesis of norovirus RNAdependent RNA polymerases by using a novel cell-based assay[END_REF]. Again no inhibition was observed at all concentrations tested (data shown in Supporting Information).

Therefore, we next hypothesized that rutin may be targeting another viral protein. Recent studies have reported that rutin is efficient for the inhibition of enterovirus 3C protease [START_REF] Lin | Fisetin and rutin as 3C protease inhibitors of enterovirus A71[END_REF], which is a picornavirus, the most closely related virus family to caliviruses. Thus, we docked rutin in the norovirus protease (PDB ID: 3UR6 [START_REF] Kim | Broad-spectrum antivirals against 3C or 3C-like proteases of picornaviruses, noroviruses, and coronaviruses[END_REF]): the predicted Kd was found to be of 588 nM (docking score of -8.5). Therefore, we performed FRET-based rutin inhibition assays with the MNV-1 and GII.4 HuNoV proteases [START_REF] Chang | Characterization and inhibition of norovirus proteases of genogroups I and II using a fluorescence resonance energy transfer assay[END_REF] at concentrations of up to 100 M. Unfortunately, no inhibition of protease activity by rutin was observed (data shown in Supporting Information).

Taken together, the in vitro assays demonstrated that the antiviral activity of rutin is not mediated by the viral polymerase or protease, suggesting the antiviral activity is mediated by one or more of the remaining viral proteins and/or unknown cellular proteins. MNV encodes nine proteins, which all could be potential targets [START_REF] Thackray | Murine noroviruses comprising a single genogroup exhibit biological diversity despite limited sequence divergence[END_REF]. However, it is unlikely that rutin targets the major capsid protein VP1, as no similarities were observed between rutin and other known small molecules targeting the VP1 capsid protein (such as suramin, heparin, citrate or fucose analogs) [START_REF] Arias | Progress towards the prevention and treatment of norovirus infections[END_REF]. Another potential explanation is that rutin is a prodrug of the active compound, for example via its sugar part (rutinose), which is synthesized in the cell. Moreover, amongst the non-nucleoside analog RdRp inhibitors are styrylchromones [START_REF] Arias | Progress towards the prevention and treatment of norovirus infections[END_REF], which with rutin a chromone moiety. Styrylchromones [START_REF] Rocha-Pereira | E)-2-styrylchromones as potential anti-norovirus agents[END_REF] are capable of inhibiting MNV replication in cell culture through an unknown mechanism, but which is hypothesized to be via targeting the RdRp. However, since the cell-based RdRp assay tests demonstrated no inhibition, it is unlikely that rutin acts as a prodrug on the MNV RdRp. At this point, we cannot rule out that rutin interacts with one of the remaining six viral proteins NS1-5, VF1 or VP2.

Time of addition experiments were performed and results are presented in Figure 3. We observed that no differences are found when 100M of rutin are added one hour before infection, at the time of infection, one hour or four hours later. However, when 1M of rutin is added, the % of virus titer drops from 76 to 39 depending if addition is made before infection or four hours after. These data lead us to hypothesize that rutin may act at the level of replication or later in the life cycle, but not at the stage of adhesion or internalization. Alternatively, the antiviral activity of rutin may be mediated through interaction with one or more cellular proteins needed during infection. This appears to be the most likely explanation as rutin inhibited MNV infection at a concentration of 1 M (Figure 2a), a concentration far below the cytotoxic concentration of rutin, 100 M (Figure 1a). However, identification of targets will require additional studies in the future.

In summary, we report that rutin, a flavonoid derivative, has antiviral efficacy against MNV-1 in cell culture. No inhibition was observed in in vitro assays with the viral polymerase or protease. The exact mechanism of action of rutin that leads to antiviral effects remains to be identified; nonetheless, we hypothesize that rutin may act at the level of replication or later in the life cycle.

Figure 1 .

 1 Figure 1. Cell viability following drug treatment. Cells were treated for 8 hours with indicated

Figure 2 . 1 (

 21 Figure 2. Antiviral efficacy of predicted RdRp inhibitors. Raw 264.7 cells were infected with MNV-1

Table 1 . Select candidates with highest docking scores a after computational screening.

 1 Docking scores are estimations of the binding free energy in kcal/mol.

	MNV	HNV	Price
	ZINC ID			Popular name
	Score b	Score c	($)/mol
	ZINC28639308 -12.1	-10.5	38	AD132644
	ZINC11726792 -9.8	-9.6	17	Sinefungin
	ZINC72206342 -9.7	-9.2	0.001	Rutin
	ZINC08143614 -9.7	-7.2	1.32	Ouabain
	ZINC04216676 -9.1	-8.4	not avail.	Kaempferitrin
	ZINC26143524 -9.1	-6.8	163	Tezosentan
	ZINC33359852 -8.9	-8.1	27	Lymecycline
	ZINC04097715 -8.8	-9.9	4	Antirrhinin
	ZINC04216682 -8.5	-10.4	0.01	Kasugamycin
	ZINC03869383 -8.4	-8.9	1	(none)
	ZINC01902827 -8.1	-9.1	15	Fidexaban
	ZINC01910616 -8.1	-8.6	5	Dabigatran
	ZINC03993831 -7.9	-8.7	81	Tomopenem
	ZINC08214483 -7.9	-10.9	0.03	Amikacin
	ZINC53684226 -7.4	-8.9	2	APAD
	ZINC04557129 -7.3	-8.2	7	(none)

a b Mean docking scores for murine norovirus (MNV) RdRp between the three chains of PDB 3UPF. c Docking scores for human norovirus (HNV) Norwalk virus RdRp from PDB 1SH2.

Acknowledgments

This work was supported by Defense Advanced Research Projects Agency (DARPA) Contract HR0011-11-C-0093. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The authors thank Audrey E. Fischer, Joshua T. Wolfe, Jeffrey S. Lin, and Andrew B. Feldman from the Applied Physics Laboratory of John Hopkins University (MD, USA) for fruitful discussions. We further thank Xiaoyan Lin and Emma Winkler at Indiana University for performing the RdRp assays and David George at Kansas State University for performing the protease assay. The computations performed for this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.