%  LaTeX support: latex@mdpi.com 
%  In case you need support, please attach all files that are necessary for compiling as well as the log file, and specify the details of your LaTeX setup (which operating system and LaTeX version / tools you are using).

%=================================================================
\documentclass[molecules,article,submit,moreauthors,pdftex]{Definitions/mdpi} 

% If you would like to post an early version of this manuscript as a preprint, you may use preprint as the journal and change 'submit' to 'accept'. The document class line would be, e.g.,https://www.overleaf.com/project/5f433acdde7a360001c5c452 \documentclass[preprints,article,accept,moreauthors,pdftex]{mdpi}. This is especially recommended for submission to arXiv, where line numbers should be removed before posting. For preprints.org, the editorial staff will make this change immediately prior to posting.

%--------------------
% Class Options:
%--------------------
%----------
% molecules
%----------
% Choose between the following MDPI journals:
% acoustics, actuators, addictions, admsci, aerospace, agriculture, agriengineering, agronomy, algorithms, animals, antibiotics, antibodies, antioxidants, applsci, arts, asc, asi, atmosphere, atoms, axioms, batteries, bdcc, behavsci , beverages, bioengineering, biology, biomedicines, biomimetics, biomolecules, biosensors, brainsci , buildings, cancers, carbon , catalysts, cells, ceramics, challenges, chemengineering, chemistry, chemosensors, children, cleantechnol, climate, clockssleep, cmd, coatings, colloids, computation, computers, condensedmatter, cosmetics, cryptography, crystals, dairy, data, dentistry, designs , diagnostics, diseases, diversity, drones, econometrics, economies, education, ejihpe, electrochem, electronics, energies, entropy, environments, epigenomes, est, fermentation, fibers, fire, fishes, fluids, foods, forecasting, forests, fractalfract, futureinternet, futurephys, galaxies, games, gastrointestdisord, gels, genealogy, genes, geohazards, geosciences, geriatrics, hazardousmatters, healthcare, heritage, highthroughput, horticulturae, humanities, hydrology, ijerph, ijfs, ijgi, ijms, ijns, ijtpp, informatics, information, infrastructures, inorganics, insects, instruments, inventions, iot, j, jcdd, jcm, jcp, jcs, jdb, jfb, jfmk, jimaging, jintelligence, jlpea, jmmp, jmse, jnt, jof, joitmc, jpm, jrfm, jsan, land, languages, laws, life, literature, logistics, lubricants, machines, magnetochemistry, make, marinedrugs, materials, mathematics, mca, medicina, medicines, medsci, membranes, metabolites, metals, microarrays, micromachines, microorganisms, minerals, modelling, molbank, molecules, mps, mti, nanomaterials, ncrna, neuroglia, nitrogen, notspecified, nutrients, ohbm, optics, particles, pathogens, pharmaceuticals, pharmaceutics, pharmacy, philosophies, photonics, physics, plants, plasma, polymers, polysaccharides, preprints , proceedings, processes, proteomes, psych, publications, quantumrep, quaternary, qubs, reactions, recycling, religions, remotesensing, reports, resources, risks, robotics, safety, sci, scipharm, sensors, separations, sexes, signals, sinusitis, smartcities, sna, societies, socsci, soilsystems, sports, standards, stats, surfaces, surgeries, sustainability, symmetry, systems, technologies, test, toxics, toxins, tropicalmed, universe, urbansci, vaccines, vehicles, vetsci, vibration, viruses, vision, water, wem, wevj

%---------
% article
%---------
% The default type of manuscript is "article", but can be replaced by: 
% abstract, addendum, article, benchmark, book, bookreview, briefreport, casereport, changes, comment, commentary, communication, conceptpaper, conferenceproceedings, correction, conferencereport, expressionofconcern, extendedabstract, meetingreport, creative, datadescriptor, discussion, editorial, essay, erratum, hypothesis, interestingimages, letter, meetingreport, newbookreceived, obituary, opinion, projectreport, reply, retraction, review, perspective, protocol, shortnote, supfile, technicalnote, viewpoint
% supfile = supplementary materials

%----------
% submit
%----------
% The class option "submit" will be changed to "accept" by the Editorial Office when the paper is accepted. This will only make changes to the frontpage (e.g., the logo of the journal will get visible), the headings, and the copyright information. Also, line numbering will be removed. Journal info and pagination for accepted papers will also be assigned by the Editorial Office.

%------------------
% moreauthors
%------------------
% If there is only one author the class option oneauthor should be used. Otherwise use the class option moreauthors.

%---------
% pdftex
%---------
% The option pdftex is for use with pdfLaTeX. If eps figures are used, remove the option pdftex and use LaTeX and dvi2pdf.

%=================================================================
\firstpage{1} 
\makeatletter 
\setcounter{page}{\@firstpage} 
\makeatother
\pubvolume{xx}
\issuenum{1}
\articlenumber{5}
\pubyear{2020}
\copyrightyear{2020}
%\externaleditor{Academic Editor: name}
\history{Received: date; Accepted: date; Published: date}
%\updates{yes} % If there is an update available, un-comment this line

%% MDPI internal command: uncomment if new journal that already uses continuous page numbers 
%\continuouspages{yes}

%------------------------------------------------------------------
% The following line should be uncommented if the LaTeX file is uploaded to arXiv.org
%\pdfoutput=1

%=================================================================
% Add packages and commands here. The following packages are loaded in our class file: fontenc, calc, indentfirst, fancyhdr, graphicx, lastpage, ifthen, lineno, float, amsmath, setspace, enumitem, mathpazo, booktabs, titlesec, etoolbox, amsthm, hyphenat, natbib, hyperref, footmisc, geometry, caption, url, mdframed, tabto, soul, multirow, microtype, tikz

%=================================================================
%% Please use the following mathematics environments: Theorem, Lemma, Corollary, Proposition, Characterization, Property, Problem, Example, ExamplesandDefinitions, Hypothesis, Remark, Definition, Notation, Assumption
%% For proofs, please use the proof environment (the amsthm package is loaded by the MDPI class).

%=================================================================
% Full title of the paper (Capitalized)
\Title{Searching Hit Potential Antimicrobials in Natural Compounds Space  Against Biofilm Formation}

% Author Orchid ID: enter ID or remove command
\newcommand{\orcidauthorA}{0000-0001-5190-4171} % Add \orcidA{} behind the author's name
\newcommand{\orcidauthorB}{0000-0003-0897-8209} % Add \orcidB{} behind the author's name
\newcommand{\orcidauthorC}{0000-0002-6492-2640}

% Authors, for the paper (add full first names)
\Author{Roberto Pestana-Nobles $^{1\ddagger}$\orcidA{}, Jorge A. Leyva-Rojas $^{1,\ddagger}$\orcidB{} and Juvenal Yosa $^{1,\ddagger}$\orcidC{}*}

% Authors, for metadata in PDF
\AuthorNames{Roberto Pestana-Nobles, Jorge A. Leyva-Rojas and Juvenal Yosa}

% Affiliations / Addresses (Add [1] after \address if there is only one affiliation.)
\address{%
$^{1}$ \quad Universidad Simón Bolivar, Facultad de Ciencias Básicas y Biomédicas, Laboratorio de Simulación Molecular y Bioinformática, Barranquilla, 080002, Colombia.}

% Contact information of the corresponding author
\corres{Correspondence: juvenal.yosa@unisimonbolivar.edu.co}

% Current address and/or shared authorship
%\firstnote{Current address: Affiliation 3} 
\secondnote{These authors contributed equally to this work.}
% The commands \thirdnote{} till \eighthnote{} are available for further notes

%\simplesumm{} % Simple summary

%\conference{} % An extended version of a conference paper

% Abstract (Do not insert blank lines, i.e. \\) 
\abstract{Biofilms are communities of microorganisms that can colonize biotic and abiotic surfaces playing a significant role in the persistence of bacterial infection and antibiotic resistance. About 65\% and 80\% of microbial and chronic infections are produced by biofilm formation. The increase in infections by multi-resistant bacteria draws attention to the discovery of new drugs based on natural inhibitory molecules. The inhibition of diguanylate cyclases (DGCs), the enzyme implicated in the synthesis of the second messenger, cyclic diguanylate (c-di-GMP), involved the biofilm formation, represents a potential method for preventing the biofilm development. It has been extensively studied using PleD protein as a model of DGC for \textit{in silico} studies as virtual screening and as a model for \textit{in vitro} studies in biofilms formation. In the present study 224205 molecules from natural products database, ZINC15  has been evaluated through molecular docking and molecular dynamic simulation, our result suggests \textit{trans}-Aconitic acid (TAA) as a possible starting point for hit-to-lead methodologies to obtain new molecules capable of inhibiting the PleD protein and hence blocking the biofilm formation.}

% Keywords
\keyword{Biofilms; virtual screening; molecular dynamics; natural products; binding energy; \textit{trans}-Aconitic acid; hit-to-lead}

% The fields PACS, MSC, and JEL may be left empty or commented out if not applicable
%\PACS{J0101}
%\MSC{}
%\JEL{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Diversity
%\LSID{\url{http://}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Applied Sciences:
%\featuredapplication{Authors are encouraged to provide a concise description of the specific application or a potential application of the work. This section is not mandatory.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Data:
%\dataset{DOI number or link to the deposited data set in cases where the data set is published or set to be published separately. If the data set is submitted and will be published as a supplement to this paper in the journal Data, this field will be filled by the editors of the journal. In this case, please make sure to submit the data set as a supplement when entering your manuscript into our manuscript editorial system.}

%\datasetlicense{license under which the data set is made available (CC0, CC-BY, CC-BY-SA, CC-BY-NC, etc.)}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Toxins
%\keycontribution{The breakthroughs or highlights of the manuscript. Authors can write one or two sentences to describe the most important part of the paper.}

%\setcounter{secnumdepth}{4}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}
Biofilms are communities of microorganisms that can colonize surfaces, they have been found on both biotic and abiotic surfaces playing a significant role in the persistence of bacterial infection. The microorganisms are embedded in a self-produced extracellular matrix of extracellular polymeric substances (EPS) \cite{Flemming2016,Yin2019}.
The biofilms formation provides advantages for the microorganisms. The change from free-living bacteria to biofilm involves the production of adhesins and extracellular matrix compound interconnecting the cells in the biofilms \cite{:Tolker-Nielsen}. The extracellular matrix helps to support and connect the cell-to-cell interaction in the biofilms and plays an important purpose in several processes including cell attachment, cell-to-cell connecting, structural function, and antibiotic tolerance. This matrix producing by the bacteria is mainly composed of proteins, enzymes, polysaccharides, signalling molecules and extracellular DNA \cite{Nadell2016}.


This microorganism assemblage provides a protected mode of growth with some resistance to antimicrobial and to killing by the host immune system. The biofilm development makes possible the surviving of the pathogen microorganisms in hostile environments and to disperse and colonize other niches also. \cite{DelPozo2018}.

In the medical field, the biofilm phenotype has been recognized as an active agent in the development of many infections, being the most common, those related to the use of medical devices\cite{fleming2017}. The National Institutes of Health have reported that, among all microbial and chronic infections, 65\%, and 80\%, respectively, are associated with biofilm formation \cite{Jamal2018}.
The understanding of the molecular and physiological mechanisms of biofilm formation will bring a tool to inhibit its establishment.

The cyclic diguanylate (c-di-GMP) is an important second messenger involved in the biofilm organization\cite{Fernicola147}. This signalling molecule modulates bacterial growth phenotypes but not only limited to biofilm formation, virulence factor production, and motility.  \cite{Fernicola147,Cai2020}. The synthesis of c-di-GMP depends on diguanylate cyclases (DGCs) which take two molecules of guanosine-5'-triphosphate (GTP) to obtain c-di-GMP in a two-step reaction. While, phosphodiesterases (PDEs) hydrolyse c-di-GMP to linear di-GMP \cite{Cai2020,ute2013}. The inhibition of DGCs as a strategy for preventing the biofilm formation represents a potential method, it has been extensively studied using PleD protein (EC: 2.7.7.65) as a model of DGC for \textit{in silico} studies, using virtual screening and as a model for  \textit{in vitro} studies in biofilms formation \cite{Fernicola147,Cai2020,ute2013,Feirer2017,pabloalviz}.

The increase in infections by multi-resistant bacteria draws attention to the discovery of new drugs. Bioactive molecules are efficient natural compounds against living organisms, they are synthesized by microorganisms, plants and animals. These substances act as a chemical defence in various competitive environments \cite{Lage2018}.
An alternative for selecting DGC natural inhibitory-molecules is the search in the database of natural products, which could supply substances, especially smalls compounds, with good binding to the target enzyme.

Since the discovery of penicillin natural products plays an important role in the discovery and development of new drugs \cite{Rossiter2017,Herrmann2017}. Between 1981 and 2014 were approved 43.6\% of anti-infective drugs and 40.7\% of anticancer agents where was based on a natural product or derivatives thereof \cite{Herrmann2017}, many natural products have become a template for drug design because of often present ligand-protein binding motifs \cite{Rodrigues2016}. Some natural products have shown antibacterial activity like the described by Nofianiet al \cite{Nofiani2020} or Emiru et al \cite{Emiru2019}. 

In the present study, an evaluation of the of diguanylate cyclase, having as a model the PleD protein from \textit{Caulobacter crescentus}  against 224205  molecules from natural resources from the ZINC15 database is carried out. Here,  virtual screening VS, molecular dynamics (MD) simulations and binding free energy calculation adopting the  Molecular Mechanics Poisson-Boltzmann Surface Area method (MM/PBSA) were used for the hit searching. A summary of the complete procedure is shown in Figure \ref{procedure_fig}. This work aimed to propose potential hit substance from natural compounds database that can be optimized to generate a promising lead candidate compounds as a competitive ligand for the diguanylate cyclase proteins.

\begin{figure}[H]
    \centering
    \includegraphics[width=15cm, height=10cm]{Images/procedure.jpeg}
    \caption{Summary of the procedure used for hit identification. An initial 20 ns of  MD was performed with the PleD protein (inactive form) in complex with GTP molecule and  Mg$^{2+}$ ions, following the crystal structure of the active form 2V0N to coupling the GTP at the catalytic site. From equilibrated trajectory starting at 1.2 ns, an average structure was produced using chimera software \cite{chimera}. A library of about 224205 natural products compounds was obtained from ZINC15 database. Bond orders, tautomeric forms, stereochemistry, hydrogen atoms and protonation states were assigned automatically by the ICM package with default parameters for each compound to prepare the ligands for docking. Ligands were docked using a VS scheme from ICM using the PleD protein average structure obtained in the above procedure. The GTP geometric space in the protein was used as a reference for the calculation of grid dimensions. 100 hits were obtained from VS procedure using the ICM energy scoring function for ranking it. From 100 hits, 35 with the best ICM energy score, were selected for 10 ns of MD simulations. For each molecule, MM/PBSA was computed using 200 frames from the last 2 ns of each MD. Finally, with the best complex (Ligand 4)-PleD, ranked from MM/PBSA binding free energy, a 110 ns of MD was done, also for GTP-PleD protein complex and PleD protein. From the above three trajectories, MM/PBSA was computed again at different windows starting from 20 ns until 60 ns, using for each window 200 frames.}
    \label{procedure_fig}
\end{figure} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and Discussion}


The diguanylate cyclase PleD protein of \textit{Caulobacter crescentus}  belongs to response regulator family, its activity is controlled by phosphorylation where two cognate kinases are involved; DivJ and PleC. PleD is required for polar differentiation in the bacterial cell cycle. Bacterial cells without functional PleD are hypermobile and fail to accomplish the swarmer-to-stalked cell transition \cite{Burton1997,Paul2004,Aldridge2003,Aldridge1999}. PleD protein contains an intrinsic nucleotide cyclase activity, which converts two molecules of GTP into 3$'$,5$'$-cyclic diguanylic acid (c-di-GMP).\cite{Paul2004}. c-di-GMP is a molecule of great interest which regulates surface-adhesion properties and motility in bacteria \cite{Jenal2004,Jenal2006}. In addition, c-di-GMP is involved in bacterial biofilm formation and persistence\cite{Tamayo2007}. Due to the fact that c-di-GMP is found exclusively in bacteria, makes it a potential target for medicinal applications. 

From the pseudo-active structure (PDB-ID: 2V0N) the mode of substrate binding as far as the position of the terminal phosphates close to the P-loop, together with two Mg$^{2+}$ ions, which are coordinated by the phosphates and two invariant carboxylates (ASP327 and GLU370). Chain A from an inactive structure (PDB-ID: 1W25), was chosen for MD and VS simulation, a GTP molecule and two Mg$^{2+}$  ions were added to the catalytic site of the protein, keeping the position of the substrate (GTP) mention above, from the pseudo-active structure. Also, ASP52 was phosphorylated, parameters for the phosphate group were obtained following the procedure mentioned elsewhere \cite{yosa2014}.


Molecular dynamics was done for the complex PleD-GTP for 20 ns. The monomer of PleD  (chain A) was chosen and GTP molecule was added to the structure in the catalytic site,  together with Mg$^{2+}$ ions and phosphorylation at ASP52. Charges and atom types, for each ligand, were assigned from the gaff force field.

The system converges at 1.2 ns of simulation(Figure S1 supplementary material). From 1.2 ns an average structure for the rest of the trajectory was computed for virtual screening (VS) simulation using Chimera package \cite{chimera} Figure \ref{pled_average_struct}. First, molecular docking was done with the GTP molecule, obtained an ICM energy score of about -27.64 kcal/mol (Table S1 supplementary material). Virtual screening was done with natural compounds library from ZINC15, having as a reference the GTP molecule from the average structure. The coordinates of the ligands were extracted from the sdf file from ZINC15 natural compounds library. Bond orders, tautomeric forms, stereochemistry, hydrogen atoms and protonation states were assigned automatically by the ICM package with default parameters \cite{Neves2012}. Atom types, charges and parameters were obtained from MMFF (Merck Molecular Force Field) force field for each ligand in the VS procedure. Geometry optimization was performed using the automatic converting procedure in ICM \cite{Neves2012}).
The ICM score for the best 100 hits molecules is summarized in table S1 supplementary material. These have ICM score energies ranging from -47.13 kcal/mol to -31.62 kcal/mol. As results from the VS procedure, the first 100 natural compounds have better ICM score energy than the GTP, Figure \ref{vs_results}. The main idea of this work is to find a compound from natural product source, which can block the active site of the PleD protein, in that way, blocking the formation of the c-di-GMP. To do that, a more accurate calculation of the binding free energy was done. MM/PBSA was used in order to compute the binding free energies for the best 35 ICM energy score molecules from VS.

\begin{figure}[H]
    \centering
    \includegraphics[width=14cm, height=10cm]{Images/average_structure.png}
    \caption{Average structure of PleD protein obtained from 20 ns of MD simulation. GTP, ASP52-phosphate moiety and Mg$^{2+}$ ions are presented in the ball and stick representation.}
    \label{pled_average_struct}
\end{figure}


\begin{figure}[H]
    \centering
    \includegraphics[width=15cm, height=9cm]{Images/vs_results.jpeg}
    \caption{Best poses for the first 100 hits ranked according to ICM energy score. In orange dot envelope representation of the space found in the average structure for GTP.}
    \label{vs_results}
\end{figure}

\subsection{Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA)}

The best 35 ligands from the VS simulations were selected to perform a 10 ns of molecular dynamics simulation. From this results, the best 6 molecules with the best binding free energy MM/PBSA were studied. For the 35 ligands, the energy values were ranked between -121.39 (ZINC04501392) and -12.84 kcal/mol (ZINC15956889), the complete list is shown in the  table S2 supplementary materials. On the other hand, the GTP molecule has a binding free energy of about -178.09 kcal/mol. The difference in energy from VS and MM/PBSA calculation lies not only in the method, but the parameter and charges used for this type of calculation (see methods).  


The six ligands with the best binding free energy according to MM/PBSA  calculations, are listed in table \ref{ligands_zinc_name} with the ZINC15 id. In Figure \ref{fig_lig_figures} a 2D structure representation for the GTP and molecules are shown. Previous studies report antibacterial activity using one of these ligands or derivatives of them, for example, the citrate can be used to functionalize and synthesize nanoparticles with antibacterial activity \cite{Khatoon2017,Choudhury2019}. In searching for new nematicidal factors from \textit{Bacillus thuringiensis} against the plant-parasitic nematode \textit{Meloidogyne incognita}, a component identified as \textit{trans}-Aconitic acid (TAA) was found. The acid showed high nematicidal activity, suggesting that TAA is specifically synthesized by the bacillus as a virulence factor, TAA acid has shown activity nematicide for \textit{Bacillus thuringiensis} bacterium \cite{nameticide}. 
 
\begin{table}[H]
    \caption{Show the ZINC15 id and the common name for the best six MM/PBSA binding free energy score compounds.}
    \centering
    \begin{tabular}{ccc}
    \toprule
    \textbf{Zinc15 id} & \textbf{Name in this work} & \textbf{Common name} \\
    \midrule
         ZINC00895081 & Ligand1 & Citrate \\
         ZINC03870145 & Ligand2 & Phosphoenolpyruvic acid \\ 
         ZINC04028701 & Ligand3 & 3-carboxy-2-(carboxymethyl)oxirane-2-carboxylate \\
         ZINC04501392 & Ligand4 & \textit{trans}-Aconitic acid \\
         ZINC19336068 & Ligand5 & bis(4-hydroxy-2-oxo-2H-chromen-3-yl)acetic acid \\
         ZINC27558828 & Ligand6 & 2,3-Bis{[(2E)-3-(4-hydroxyphenyl)-2-propenoyl]oxy}succinic acid \\
    \bottomrule
    \end{tabular}
    \label{ligands_zinc_name}
\end{table}

\begin{figure}[H]
    \centering
    \includegraphics[width=12cm, height=18cm]{Images/ligandos_gtp.png}
    \caption{2D structure of the six ligands with the best MM/PBSA performance.}
    \label{fig_lig_figures}
\end{figure}  

The MM/PBSA calculation brings the information of how strong is the binding between protein and ligand, for this, the algorithm considers the binding free energy equal to the free energy of the complex minus the sum of the free energy of the protein plus the free energy of the ligand, see equation (\ref{eq_mmpbsa_general}). 

\begin{equation} \label{eq_mmpbsa_general}
\Delta G=G_{complex}-(G_{protein}+G_{ligand})
\end{equation}

Each $G$ component from equation (\ref{eq_mmpbsa_general}) can be calculated as:

\begin{equation} \label{eq_g_component}
\Delta G=\langle E_{MM}\rangle -TS+\langle G_{solvation}\rangle
\end{equation}

Where $\langle E_{MM}\rangle$ is the average molecular mechanics potential energy in vacuum, T is the temperature, S is the entropy and $\langle G_{solvation}\rangle$ is the free energy of solvation. The term $\langle E_{MM}\rangle$ include energy of both bonded and nonbonded interactions, $E_{bonded}$ and $E_{nonbonded}$. The $E_{bonded}$ are consider of bond, angle, dihedral and improper interactions. $E_{bonded}$ consider electrostatic ($E_{elec}$) and van der Waals ($E_{vdW}$) interactions, see equation (\ref{eq_emm}). $\langle G_{solvation}\rangle$ are included polar and non polar free energies, polar contribution is calculated using Poisson-Boltzmann  equation  and nonpolar contribution includes repulsive and attractive forces between solute and solvent generated by cavity, $E_{cavity}$ formation and van der Waals interactions, $E{_vdW}$, see equation (\ref{eq_sol}), for more information about how the equation works and how it is calculated consult the work of Kumari et al. \cite{Kumari2014}, Baker et al. \cite{Baker2001}, Wisser et al. \cite{weisser1999} and Konecny et al. \cite{konecny2012}.

\begin{equation} \label{eq_emm}
\langle E_{MM}\rangle = E_{bonded} + E_{nonbonded} = E_{bonded}+E_{vdW}+E_{elec}
\end{equation}

\begin{equation} \label{eq_sol}
\langle G_{solvation}\rangle = G_{polar}+G_{nonpolar} = G_{polar} + E_{cavity} + E_{vdW}
\end{equation}

A lower value of $\Delta$G means a better coupling between protein and ligand. The table \ref{table_mmpbsa} shows the values obtained for the 6 ligands with the best binding free energy. Results show that non of the molecules tested get the same or close binding free energy as GTP, the difference between GTP and ligands energies are; Ligand1 78 kcal/mol, Ligand2 74.77 kcal/mol,  Ligand3 58.51 kcal/mol, Ligand4 56.76 kcal/mol, Ligand5 70.16 kcal/mol and Ligand6 77.11 kcal/mol, the biggest difference is presented by Ligand1 and the lowest difference is with Ligand4. In order to have a better insight of which residues contribute more to the $\Delta G$ binding free energy, energy decomposition was performed (see methodology section for details), with this, the total $\Delta G$ energy is presented as a contribution per residue, with this information we can classify residues into favourable and non-favourable interactions. 

\begin{table}[H]
    \caption{Contribution of energies in the MM/PBSA calculation, all energy values are in kcal/mol.}
    \centering
    \begin{tabular}{cccccc}
    \toprule
    \textbf{Ligand} & \textbf{VDWAALS} & \textbf{EEL} & \textbf{EPB} & \textbf{ENPOLAR} & \textbf{DELTA G}\\
    \midrule
        Ligand1 & 2.92 $\pm$ 4.30  & -867.68 $\pm$ 21.67 & 766.94 $\pm$ 18.20 & -1.64 $\pm$ 0.06 & \cellcolor{lightgray}-99.46 $\pm$ 5.67 \\
         Ligand2 & 19.03 $\pm$ 4.87 & -736.68 $\pm$ 24.60 &  615.55 $\pm$ 21.83 & -1.22 $\pm$ 0.06 & \cellcolor{lightgray}-103.32 $\pm$ 4.81 \\
         Ligand3 & 8.72 $\pm$ 5.14 & -776.77 $\pm$ 38.28 & 650.15 $\pm$ 35.30 & -1.67 $\pm$ 0.06 & \cellcolor{lightgray}-119.58 $\pm$ 7.11 \\
         Ligand4 & 7.86 $\pm$ 4.83 & -807.77 $\pm$ 33.78 & 680.01 $\pm$ 29.12 & -1.43 $\pm$ 0.09 & \cellcolor{lightgray}-121.33 $\pm$ 9.13 \\
         Ligand5 & -2.20 $\pm$ 5.73 & -686.56 $\pm$ 28.76 & 583.50 $\pm$ 27.10 & -2.67 $\pm$ 0.13 & \cellcolor{lightgray}-107.93 $\pm$ 8.92 \\
         Ligand6 &  -20.27 $\pm$ 4.39 & -533.33 $\pm$ 19.94 & 456.46 $\pm$ 16.33 & -3.84 $\pm$ 0.11 & \cellcolor{lightgray}-100.98 $\pm$ 7.20 \\
        GTP & -10.86 $\pm$ 6.45 & -1287.03 $\pm$ 35.53 & 1123.46 $\pm$ 29.55 & 3.49 $\pm$ 0.17 & \cellcolor{orange}-178.09 $\pm$ 10.82 \\
    \bottomrule
    \end{tabular}
    \label{table_mmpbsa}
\end{table}

\begin{figure}[H]
    \centering
    \includegraphics[width=15cm, height=17cm]{Images/mmpbsa_gtp_lig1-2.png}
    \caption{Interaction per residue, from residue 300 to 457, all energies were calculated in kcal/mol. \textbf{(a)} Contribution of all residues on the MM/PBSA energy for GTP. \textbf{(b)} Zoom from residue 300 to 457 for the GTP, \textbf{(c)} Ligand1, \textbf{(d)} Ligand2. Residues with a contribution higher than 1 or lower than -1 are marked.}
    \label{fig_mmpbsa_gtp_lig1-2}
\end{figure} 


\begin{figure}[H]
    \centering
    \includegraphics[width=15cm, height=17cm]{Images/mmpbsa_lig3-6.png}
    \caption{Interaction per residue from residue 300 to 457, all energies were calculated in kcal/mol. \textbf{(a)} Ligand3, \textbf{(b)} Ligand4, \textbf{(c)} Ligand5 and \textbf{(d)} from Ligand6. Residues with a contribution higher than 1 or lower than -1 are marked.}
    \label{fig_mmpbsa_lig-3-6}
\end{figure}

\begin{figure}[H]
    \centering
    \includegraphics[width=15cm, height=20cm]{Images/ligands_sites_all.png}
    \caption{Contribution per residues around 3 {\AA} from the ligand in the MM/PBSA calculation for each ligand and GTP, blue colour means more favourable interaction meanwhile red colour means non-favourable interaction, the ligand is highlighted in green.}
    \label{fig_mmpbsa_site_all}
\end{figure}


Figure \ref{fig_mmpbsa_gtp_lig1-2} only take into account the residues between 300 and 457, because as is shown in Figure \ref{fig_mmpbsa_gtp_lig1-2}a the first 300 residues do not show contribution to the binding free energy. For Figure \ref{fig_mmpbsa_gtp_lig1-2} and Figure \ref{fig_mmpbsa_lig-3-6} negative values are considered as favourable interaction, the interactions between two of the three magnesium present in the protein play a role coordinating the carboxyl groups in all 6 ligands including GTP which help to gain more stability. Taking the GTP as a reference, an interaction is observed with ARG445, this is a highly favourable interaction which is not present in any of the 6 ligands evaluated, suggesting this residue as an important residue in the binding pocket which help to reach a more stable conformational configuration. The ligands show common interactions with LYS331, LYS441, PHE329 and PHE330 as favourable except for the Ligand2 presumably because this ligand is the shortest one limiting its interaction with the protein.  
The result of positive energy suggest a non-favourable interaction with the protein residues, all ligands show non-favourable interaction with ASP326, ASP328, GLU369, ILE327 and GLU370, these four residues are common and have the higher positive values among the ligands. All ligands except for Ligand4 and Ligand6 show non-favourable interaction with ASP434. The non-favourable interaction for the GTP is mainly with ASP326 and GLU369.


\subsubsection{Root Mean Square Deviation (RMSD)}

From the 10 ns molecular dynamic simulation used for the calculation of MM/PBSA, a RMSD analysis was done. The RMSD was used to verify the protein stability during the entire simulation \cite{rmsd}, less deviation in these value means more stability in the protein. For these, the table \ref{table_lig_rmsd} show the average of the RMSD for the 6 ligands and the GTP. As the Figure \ref{fig_rmsd_compare} show the 6 ligands stabilize around 7 ns of simulation and GTP around 4 ns. These values were obtained using cpptraj tool \cite{cpptraj} aligning to the first frame of the production part of the simulation, the result suggests a stable configuration of the protein with the ligands as we can see the same behaviour in the Figure \ref{fig_rmsd_compare}, GTP shows less variation compared to the rest of the ligands. 

\begin{table}[H]
    \caption{The average RMSD for each Pled-ligand and Pled-GTP complexes, in angstroms with its standard deviation considering the entire production.}
    \centering
    \begin{tabular}{ccc}
    \toprule
    \textbf{Ligand} & \textbf{RMSD} &\textbf{Difference between GTP and ligand}\\
    \midrule
         Ligand1 & 2.261 $\pm$ 0.365 & 0.284\\
         Ligand2 & 2.189 $\pm$ 0.336 & 0.212\\
         Ligand3 & 2.165 $\pm$ 0.411 & 0.188\\
         Ligand4 & 2.480 $\pm$ 0.255 & 0.503\\
         Ligand5 & 2.506 $\pm$  0.350 & 0.529\\
         Ligand6 &  2.305 $\pm$  0.389 & 0.328\\
         \rowcolor{lightgray}
        GTP & 1.977 $\pm$ 0.289 & 0\\
    \bottomrule
    \end{tabular}
    \label{table_lig_rmsd}
\end{table}

\begin{figure}[H]
    \centering
    \includegraphics[width=15cm, height=5cm]{Images/rmsd_all_ligands_gtp.png}
    \caption{RMSD calculation for the 10 ns of production for each ligand and GTP.}
    \label{fig_rmsd_compare}
\end{figure}  

\subsubsection{Hydrogen bonds (H-bonds)}

Hydrogen bonds represent an important interaction between protein-ligand complexes and the main reason for protein-ligand selectivity, this is why hydrogen bond analysis is often used to describe affinity between protein and ligand \cite{Nittinger2017}. 
Using hbond function in the cpptraj tool we analyze the average amount of time at which hydrogen bond was present during the simulation, from this information, it can be seen, which ones among all the residues, in the protein have more H-bond formation with the ligands. For these, the H-bond analysis for each ligand (including GTP) was done using a threshold of 3 {\AA} around the ligands during all the production taking in account frame by frame. 

\begin{table}[H]
    \caption{Residues involved in the H-bond formation with the percentage of how long this bond was present during the entire 10 ns of simulation.}
    \centering
    \begin{tabular}{cc}
    \toprule
    \textbf{Ligand} & \textbf{Residues involve in the h-bond formation with the major contribution}\\
    \midrule
        Ligand1 & LYS441(54.20\%), LYS331(35.88\%), PHE330(7.11\%), PHE329(2.80\%) \\
         Ligand2 & LYS441(85.28\%), PHE329(9.37\%), PHE330 (4.97\%) \\
         Ligand3 & LYS441(38.70\%), PHE329(24.17\%), PHE330(19.74\%), LYS331(17.39\%)\\
         Ligand4 & LYS441(52.68\%), LYS331(24.05\%), PHE330(19.28\%), PHE329(3.99\%) \\
         Ligand5 & LYS441(54.31\%), PHE330(35.45\%), LYS331(9.31\%) \\
         Ligand6 &  LYS441(36.13\%), LYS331(23.48\%), PHE330(17.81\%), ASN334(17.03\%), LYS332(4.77\%) \\
         \rowcolor{lightgray}
        GTP & ARG445(39.05\%), LYS331(34.63\%), LYS441(17.35\%), PHE330(8.12\%) \\
    \bottomrule
    \end{tabular}
    \label{table_hbond}
\end{table}


Among all ligands evaluated, the LYS441 is the more important residue in the H-bond formation. In some cases being present in more than the 80\% of the production (Ligand2). Residues PHE331, PHE330 and LYS331 were also present as relevant. For the GTP the residue with more interactions was the ARG445, this also is seen the MM/PBSA decomposition results.

A second analysis was performed to observe the behaviour of the best ligand obtained in the MM/PBSA calculation, which is Ligand4. For this, 100 ns of molecular dynamics with the same condition of the previous 10 ns were done. The starting point for this new simulation was the ending of the previous one giving a total of 110 ns of simulation. In this second analysis, the PleD-protein without any ligand, to serve as a reference, the complex PleD-GTP and the complex PleD-Ligand4, were taken for simulations.  MM/PBSA calculations were done, using different time windows,  20, 30, 40, 50 and 60 ns of the molecular simulation were consider, in each window, the same amount of snapshots (200) were taken.


\subsection{\textit{Trans}-Aconitic Acid (Ligand4)} 

\begin{figure}[H]
    \centering
    \includegraphics[width=14cm,height=20cm]{Images/second_analysis_hbond_rmsd_rmsf_3.png}
    \caption{Analysis results  from the second production of 110 ns, \textbf{(a)} RMSD in angstrom, \textbf{(b)} RMSF in angstrom, \textbf{(c)} Radius of gyration (Rg) in angstrom, \textbf{(d)} Hydrogen bonds for GTP and \textbf{(e)} Hydrogen bonds for Ligand4.}
    \label{second_analysis_hbond_rmsd_rmsf}
\end{figure}


\subsubsection{Root Mean Square Deviation (RMSD)}

The RMSD show high stability in the three cases (PleD, GTP and Ligand4), \ref{second_analysis_hbond_rmsd_rmsf}.a, the average value with it standard deviation for the PleD-protein 2.75 $\pm$ 0.31, PleD-GTP complex 2.57 $\pm$ 0.28 and PleD-Ligand4 complex 3.20 $\pm$ 0.45. The PleD-protein present a higher value than PleD-GTP meaning GTP help the protein itself to reach a more stable configuration. The PleD-Ligand4 is the one who presents more variation although has high stability during the dynamic as is represented in the low standard deviation.

\subsubsection{Root Mean Square Fluctuation (RMSF)}

For the RMSF the same approach as the MM/PBSA was made, the graph just take into account residues 300-457 which are the residues were the interaction are establish, Figure \ref{second_analysis_hbond_rmsd_rmsf}b, the same fluctuation is observed in the protein structure for the three cases, specifically in the PleD-GTP, a major variation in residues 450-455 is present, the average variation with the standard deviation is 1.65 $\pm$ 0.63 for the PleD-protein, 1.53 $\pm$ 0.87 for PleD-GTP complex and 1.53 $\pm$ 0.62 for PleD-Ligand4 complex. The RMSF analysis agrees with the RMSD showing the PleD by itself present lower stability compared to a binding PleD with GTP.

\subsubsection{Radius of gyration}

The Radius of gyration (Rg) can give an idea of the compactness of a protein structure \cite{Lobanov2008}, the Figure \ref{second_analysis_hbond_rmsd_rmsf}c suggest structural stability in the 3 systems here considered however the complex PleD-Ligand4 less Rg value in comparison of the other two structures. The average value with the standard deviation for the three systems are 25.48 $\pm$ 0.27 for PleD-protein, 25.36 $\pm$ 0.26 PleD-GTP and 24.1 $\pm$ 0.17 PleD-Ligand4. 



\subsubsection{Hydrogen bond (H-bond)}

The H-bond analysis results for the 110 ns production (\ref{table_hbond_110ns}) show differences compared to the results from the table \ref{table_hbond}, for the GTP is observed the same tendency as the previous result, where the residue ARG445 is the one which is present for more time during the simulation, in this case, the order is maintained, the LYS331, LYS441 and PHE330 residues come after the ARG445. In the case of the Ligand4 there are different order, in this case, the residue PHE330 is the one with more presence followed by LYS331 and LYS441. 


\begin{table}[H]
    \caption{Residues involved in the H-bond formation with the percentage of how long this bond was presented during the entire simulation for the 110 ns of simulation.}
    \centering
    \begin{tabular}{cc}
    \toprule
    \textbf{Ligand} & \textbf{Residues involve in the h-bond formation with the major contribution}\\
    \midrule
        PleD-Ligand4 & PHE330(9.7\%) ,LYS331(6.01\%),  LYS441(5.69\%) \\
        PleD-GTP & ARG445(41.44\%), LYS331(15.14\%), LYS441(12.07\%), PHE330(11.95\%)  \\
    \bottomrule
    \end{tabular}
    \label{table_hbond_110ns}
\end{table}

The H-bond formation shows that, during the production, the GTP is making even nine H-bonds while the Ligand4 just make only four as is shown in Figure \ref{second_analysis_hbond_rmsd_rmsf}d and \ref{second_analysis_hbond_rmsd_rmsf}e. The GTP present the double of H-bond compared to the Ligand4, this is due to the phosphate group (the GTP has 3 of them) which provide more possible points to establish H-bonds thanks to the hydroxyl group (10 in the whole structure), instead of the Ligand4 which only have 6 hydroxyl group.  

\subsubsection{Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA)}

The MM/PBSA values here presented are averages of the windows establish in this step, the values are (with the standard deviation); -169.93 $\pm$ 4.03 for PleD-GTP and -98.79 $\pm$ 0.43 for PleD-Ligand4. In the Figure \ref{second_analysis_mmpbsa} are represented the average values of the residues from the different windows (20, 30, 40, 50 and 60), compared to the previous results (Figure \ref{fig_mmpbsa_gtp_lig1-2}) for the GTP, here can be seen a similarity, which indicates the same tendency for the residues like ASP326, ASP328 and GLU369 which have non-favourable interaction while ARG445, LYS441, LYS331, PHE330 and PHE329 have the major favourable contribution to the interaction with the GTP, suggesting the importance of these residues in the binding pocket. In the case of Ligand4, the residues showed the same tendency as in the simulation of 10 ns, ASP326, ILE327, ASP328, GLU369 and GLU370 residues, show non-favourable interaction while residues PHE329, LYS331 and LYS441 show favourable interaction. In both cases, the magnesium atoms present an important contribution to the binding free energy. 

\begin{figure}[H]
    \centering
    \includegraphics[width=13cm]{Images/second_analysis_mmpbsa.png}
    \caption{MM/PBSA Average ligand binding free energies for the 110 ns of simulation. Here, the residues 300 to 457 were taken into account. \textbf{(a)} GTP and \textbf{(b)} Ligand4.}
    \label{second_analysis_mmpbsa}
\end{figure}

The low standard deviation for the reported result at different windows is interpreted as high stability of the complexes, comparing the two MM/PBSA result (first and second analysis) some variation in energy are presented and the emergence or disappearance of new residues as GLU370 or ASP343 for GTP, and LYS332 for the Ligand4, is due to better conformational sampling obtained with more simulation time (110 ns). The result of the MM/PBSA with a short dynamic gives us an idea of the overall interaction and also of possible favourable or non-favourable ones. 

For the GTP results, where the major contribution for the MM/PBSA binding free energy is the H-bond formation from ARG445,  suggest this residue as a key one for the ligand-protein interaction. The absence of this interaction with the rest of the ligands can indicate a lower binding free energy in comparison with the GTP. Both MM/PBSA results, the one obtained by 10 ns and the average from the five windows, showed that the Ligand4 has good stability, suggesting this molecule as a possible template for design new ligand who can surpass the GTP binding free energy and effectively inhibit the PleD protein. 

In this research the binding free energy was calculated by MM/PBSA method and was used to define a possible compound capable of inhibiting the PleD protein and hence blocking the formation of biofilms, the \textit{trans}-aconitic acid performed the best among the 224205 ligands here studied. This compounds is found in plants like \textit{Zea mays}, \textit{Triticum aestivum}, \textit{Avena sativa}, \textit{Brachiaria plantaginea} and \textit{Saccharum officinarum} \cite{Bortolo2018}. Methodologies of extraction are described by Schnitzler et al. \cite{Schnitzler2007} and Kanitkar et al. \cite{kanitar2013}. In our results \textit{trans}-aconitic acid does not surpass GTP binding free energy, however it does not mean that can not be used, our results suggest good interaction with the PleD protein, further optimization can be made to improve the results, strategies like Fragment-Based Drug Discovery (FBDD) \cite{deSouzaNeto2020} or Fragment Molecular Orbital (FMO) \cite{KITAURA1999701, HEVENER2018265} can be use to optimized the hit to become an effective lead.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Materials and Methods}

\subsection{Data collection}

To search for antimicrobial agents from natural compounds, the ZINC15 \cite{Sterling2015} database was chosen. ZINC15 is a database with public access established to allow ready access to compounds for virtual screening. Molecules were downloaded in SDF format. ZINC 15 has a natural product-related compounds database of about 224205. All records from naturals products were download and stores in a single file.

\subsection{Virtual screening (VS)}

Structure-based virtual screening applying docking simulations was performed using ICM software package \cite{abagyan_ICM,abagyab_ICM_2,toreov_ICM,Jianghong_ICM,Fernandez_ICM,Fernandez_ICM_2,mendez_ICM,mendez_ICM2,Neves2012}. Virtual Ligand Screening (VLS) in ICM was performed by docking the database downloaded from ZINC15 against PleD protein structure obtained from protein data bank PDB (PDB:ID 1W25) followed by an evaluation of the docked conformation with a binding-score function.  The search space encompassed the catalytic domain, the DGC (GGDEF) domain consists of a five-stranded central $\beta$-sheet surrounded by $\alpha$ helices (Figure \ref{pled_average_struct}). Water molecules and ligands had been removed from the original crystal structure. Phosphoryl analog BeF$^{-3}$ attached to ASP52, was changed by phosphate group which was parameterized using the same procedure as was mention elsewhere \cite{yosa2014}. A grid for space search was set up using as a reference the GTP$\alpha$S and Mg$^{2+}$ from 2V0N crystal structure, with the following dimensions in {\AA}: centre (x, y, z) = (52.14.6, 76.2), dimensions (x, y, z) = (21.2, 17.9, 17.9). The docking simulation was then run at exhaustiveness of 5 and set to only output the lowest energy pose.

\subsection{Ligand binding free energy calculations}

For the ligands, no geometry optimization was applied because we kept the orientation obtained from the docking results. Charge for ligands were calculated with \textit{ab initio} methods using B3LYP and 6-31G as basis set in GAUSSIAN09 software \cite{gausssian09}. With the output file from gaussian09, charges and topology file (prepi and frcmod AMBER files respectively) were prepared using antechamber and parmchk2 from Ambertools18 \cite{amber18}. Although antechamber generates a file with the charge information of ligands these charges were changes for the ones calculated with gaussian. 

\subsection{System Preparation}

The system consists of the  protein-ligand complex embedded in a rectangular box with waters and ions,  LEaP module of Amber was used for the preparation of all  systems \cite{amber18}. Initially, the mbondi2 parameter for PBRadii were set \cite{Su2015}. The ff14SB force field was used for the protein which is an improvement of ff99SB \cite{Simmerling2015}. 35 ligands with the best ICM energy score were chosen. The GAFF force field was used for the ligands\cite{Wang2004}, all the system were embedded in a box of size 100, 79, 109 {\AA} (approximately)  for X, Y, Z respectively, we use TIP3P water type molecules \cite{water_new,water_old} and Na$^{+}$ ions were add to neutralize the charge. 

\subsection{Molecular Dynamics (MD)}

For all systems, first, a minimization of only waters molecules was done, for this, we restrain all the solute (protein and ligand) with a force of 100 kcal/mol-\r{A}$^2$, in this minimization we establish the maxcyc and the ncyc to 50000 and 1500 respectively, the ntc was 1 (the shake was off) this is recommended to do it during the minimization \cite{amber18}, constant volume was established. Then a minimization of the entire system using this time maxcyc=100000 and ncyc=1000, again constant volume was done. 

After energy minimization the system was gradually heated until 300 K for 50 ps, in this part the solute was restrained with a force of 2.0 kcal/mol-\r{A}$^2$, for the last step the system was equilibrated for 5 nanoseconds (ns), we establish ntt=3, ntp=1, ntc=2 and ntf=2, with constant pressure. Finally, for the first analysis, we perform the production for 10 ns in NPT ensemble. 

The second analysis of PleD protein, PleD-GTP and PleD-Ligand4 complexes, a 100 ns more were simulated. For the PleD-GTP and PleD-Ligand4, the 100 ns production start at the end of the previous 10 ns, using the same conditions as the 10 ns production.   

\subsection{Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) Calculation}

The binding free energy was calculated using the MM/PBSA script from AMBER package \cite{Miller2012}, in the first analysis a total of 200 snapshots were taken for the MM/PBSA analysis, these snapshots are equal to the last 2 ns of simulation, Su et al suggest that working with data of MD trajectories longer than 2 ns decreases the accuracy of the prediction \cite{Su2015}. The parameters igb = 2, inp = 1 and radiopt = 0 were used, the igb parameter is chosen depending on the mbondi parameters following the amber manual indications, the inp uses the solvent-accessible-surface area to correlate total nonpolar solvation free energy and radiopt indicating to use the default atomic radius of the input file \cite{amber18}. In the calculation of MM/PBSA binding free energy entropy was not considerate, this due to the high computational cost of it and this contribution just improve the result slightly or even can decrease the effectiveness of the MM/PBSA method \cite{Ben-shalom,gen2011,Hou2011}. For the decomposition per residue of the MM/PBSA result, the \&decomp general option was established with idecomp=1 in the same input file.  

In the second analysis, we take the last 20, 30, 40, 50, 60 ns of production and  200 snapshots,  with the same condition as the first MM/PBSA calculation.  

\subsection{RMSD, RMSF, Radius of gyration and H-bond}

The RMSD (Root Mean Square Deviation), the RMSF (Root Mean Square Fluctuation), the radius of gyration (Rg) and H-bond (hydrogen bond formation) were calculated with cpptraj tool \cite{cpptraj}, we use the entire production to make this analysis and the first frame was taken as reference for alignment in the RMSD and RMSF analysis.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}

Computational methods, using virtual screening, play a significant role in drug discovery. The identification of the drug interaction with the target macromolecule can help to clarify the potential molecule$'$s action mechanisms.

In our computational study, the binding free energy between PleD protein and ligands was calculated through MM/PBSA to rank the best molecules obtained from virtual screening of naturals products from ZINC15 database, we identify ARG445 as a key residue in the binding pocket, residues LYS441, LYS331, PHE 330 and PHE329 present favourable interaction and GLU369 present a common non-favourable interaction among the ligands here studied. Our results suggest \textit{trans}-aconitic acid as a possible starting point for hit-to-lead methodologies. 





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{6pt} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% optional
\supplementary{The following are available online at \linksupplementary{s1}, Figure S1: title, Table S1: title, Video S1: title.}

% Only for the journal Methods and Protocols:
% If you wish to submit a video article, please do so with any other supplementary material.
% \supplementary{The following are available at \linksupplementary{s1}, Figure S1: title, Table S1: title, Video S1: title. A supporting video article is available at doi: link.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\authorcontributions{Conceptualization, R.P.N, J.A.L.R and J.Y.; methodology, R.P.N and J.Y.; software, R.P.N and J.Y.; formal analysis, R.P.N and J.Y.; writing--original draft preparation, R.P.N.; writing--review and editing, R.P.N, J.A.L.R and J.Y.; visualization, R.P.N.; supervision, J.A.L.R and J.Y.; }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\funding{This work was supported by "Convocatoria del Fondo de Ciencia, Tecnología e Innovación del Sistema General de Regalías para la conformación de una lista de proyectos elegibles para ser viabilizados, priorizados y aprobados por el OCAD en el marco del Programa de Becas de Excelencia Doctoral del Bicentenario 2019" which make possible the PhD studies of Roberto Pestana-Nobles.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgments{This work was supported by "Convocatoria del Fondo de Ciencia, Tecnología e Innovación del Sistema General de Regalías para la conformación de una lista de proyectos elegibles para ser viabilizados, priorizados y aprobados por el OCAD en el marco del Programa de Becas de Excelencia Doctoral del Bicentenario 2019" which make possible the PhD studies of Roberto Pestana-Nobles.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\conflictsofinterest{The authors declare no conflict of interest.} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\abbreviations{The following abbreviations are used in this manuscript:\\

\noindent 
\begin{tabular}{@{}ll}
EPS & Extracellular polymeric substances \\
c-di-GMP & Cyclic diguanylate \\
DGC & Diguanylate cyclases \\
PDEs & Phosphodiesterases \\
GTP & Guanosine-5'-triphosphate \\
TTA & \textit{trans}-aconitic acid\\
VS & Virtual screening \\
MD & molecular dynamics \\
MM/PBSA & Molecular Mechanics Poisson-Boltzmann Surface Area \\
MMFF & Merck Molecular Force Field \\

\end{tabular}}


% Citations and References in Supplementary files are permitted provided that they also appear in the reference list here. 

%=====================================
% References, variant A: internal bibliography
%=====================================
\reftitle{References}
\externalbibliography{yes}
\bibliography{bio.bib}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}

