CONFORMATIONAL STUDY OF MOLECULES IN A BIOLOGICAL ENVIRONMENT, DESIGN OF INHIBITORS OF HUMAN AMINOPEPTIDASE M1 IMPLICATED IN CANCER THERAPY
Issouf Soro1 , Hermann N’Guessan1 , Akoun Abou2 , Raymond Kre N’Guessan1 ,
1Laboratory of Fundamental and Applied Physics, University of Abobo Adjamé (Now Nangui Abrogoua), Côte d’Ivoire.
2Department of Training and Research in Electrical and Electronic Engineering, Research Team: Instrumentation, Image and Spectroscopy, National Polytechnic Institute, Yamoussoukro, Côte d’Ivoire.
3Laboratory of Structural and Theoretical Organic Chemistry, University of Cocody (Now Félix Houphouët Boigny), Côte d’Ivoire.
4Laboratory of Material Sciences, the Environment and Solar Energy, University Félix Houphouet Boigny, Ivory Coast.
5Quantitative Life Science, ICTP-UNESCO, Strada Costiera 11, I 34151 Trieste, Italy.
Objective: A novel subnanomolar anticancer hydroxamic acid containing drug candidates, inhibitors of human M1 aminopeptidase (APN) a recent validated target and has reached the predicted subnanomolar range of inhibitory potency.
Methods: A quantitative structure activity relationships (QSAR) complexation model has been developed from a compounds of 37 hydroxamic acid derivatives (AHD1-37 as training set, TS) to establish a linear correlation between the calculated relative Gibbs free energies (GFE: ΔΔGcom) of APN-AHDx complex formation and the experimental inhibition potency (Kiexp). The predictive power of the QSAR model was then validated first with 9 other AHDs not included in the TS and thereafter with the generation of a 3D-QSAR-PH4 pharmacophore (PH4) model to screen the AHD chemical subspace built as a virtual combinatorial library of more than 58,644 AHD analogs). Finally the best PH4 hits were evaluated with the initial QSAR model for predicted potency (Kipre) and pharmacokinetic profile.
Results: The QSAR model linear correlation equation: pKiexp=-0.1901×∆∆Gcom + 8.2886, R2=0.94, the subsequent PH4 model linear correlation between experiment and PH4-estimated Ki: pKiexp=1.0006× pKipre + 0.0028, R2=0.79 documents the high predictive power of this approach. Finally the screening of the virtual library of AHD analogs yielded 95 orally bioavailable candidates the best reaching a predicted potency (Kipre) of 50 pM and displaying favorable pharmacokinetic profile.
Conclusion: The combined use of molecular modeling (QSAR) and in silico PH4-based screening of the hypothetical combinatorial library has resulted in proposed and predicted potent anticancer candidates with a suitable pharmacokinetic profile.
Keywords: ADMET, complexation model, Drug design, molecular modeling, pharmacophore model, QSAR model.
INTRODUCTION
Cancer is one of the most worrying public health concerns in the world today. According to Global Cancer (GLOBOCAN) 2020 studies, more than 19.3 million new cases and 10 million cancer-related deaths were estimated1. Although cancer survival rates are expected to improve and cancer mortality rates have declined, cancer remains a leading cause of death worldwide. The undesirable side effects of many cancer drugs mainly are due to low selectivity towards non-cancerous cells2 and long-term use, inevitably is accompanied by drug resistance and reduction of their efficiency3. Then extensive research has been conducted to identify and characterize diverse cancer therapeutic targets at the molecular level4.Aminopeptidase N (APN/CD13) is one of the most studied cancer therapeutic targets5,6; an enzyme omnipresent in human body with multipurpose enzymatic functions, receptor for others human viruses (e.g. coronaviruses)7. Thus, APN is inescapable in the regulation of protein turnover in almost all of the organisms8-12, it dysregulation occurs in practically all types of human malignancies13. APN has been observed in several types of cancer14-18. APN activity impact metastasis, a complex biological process that concourse to more than 90% of cancer-related deaths19,20. APN is a Metallo aminopeptidase in the family of M121, characterized by a single zinc ion, involved in zinc binding and a motif substrate recognition GXMEN21,22,23,24. The x-ray crystal structure of human APN was reported by Wong et al.,21,25.Among the APN inhibitors anticancer candidates, bestatin (1) is the most studied competitive one26 against gastric, lung cancer and acute myeloid leukemia, and acute non-lymphatic leukemia27-33. Another APN inhibitor, Tosedostat (2) is an orally bioavailable prodrug converted to a pharmacologically active drug inside cells34-38. Later Jisook Lee et al., repurposed the compound (3) N-(2-(hydroxyamino)-2-oxo-1-[3',4',5'-trifluoro(1,1'-biphenyl)-4-yl]ethyl)pivalamide reporting it as a novel APN inhibitor more potent than bestatin (1) and Tosedostat (2). Moreover they synthesized a series of hydroxamic acid inhibitors (4) to optimize binding interactions around and beyond the S1' subsite of APN, the most potent being compound (5) in the low nanomolar range, Kiexp=4.5 nM39.
In this work, QSAR ‘complexation’ model was built, starting from crystal structure of APN-AHD1 complex (Ki=4.5 nM, PDB entry 4FYR)39. In active site the key ligand-receptor interactions of APN-AHD1 complex shown on Figure 2 in 2D scheme were considered. Gibbs free energies of ligand-receptor complex formation (∆∆Gcom) were calculated for the series of molecules and correlated them with the observed biological activities. The resulting quantitative structure-activity relationships model (QSAR), which employs the computed parameter ∆∆Gcom was able to explain approximately 94% of the variation in the observed Kiexp. The QSAR model allowed structure-based design of novel AHD analogs. The identified virtual hits reached predicted inhibitory activities Kipre against the APN in the sub nanomolar concentration range. Metrics describing interactions at the active site of APN were assessed from analysis of the X-ray crystal structure of APN (PDB code 4FYR) in complex with one of the most active inhibitors studied in this work (5)39. The catalytic zinc binding group in the active site is coordinated by a catalytic triad His388, His392 and Glu411 (not shown in the 2D diagram in Figure 2), and the S1 pocket with Asn350, Ala351, Arg363, Gln857 Asp858, Thr860, Ser861, Phe896 and Ser897. Also, a deep hydrophobic pocket S1' with residues Arg381, Ser415, Glu419, Tyr419 Arg442 and Tyr477.
METHODS
Training and validation sets
The literature had been used for training and validation sets inhibitors of hydroxamic acid analogs of human APN39. Their Kiexp covers a very wide range (4.5 ≤ Kipre ≤ 4,420 nM), more than four orders of magnitude, suitable for a reliable QSAR model. Out of a total of 46 compounds, 37 were used for the training set (TS) and 9 for the validation set (VS).
Model building
Three-dimensional (3D) molecular models of free inhibitors (I), free APN enzyme (E) and enzyme-inhibitor complexes (E:I), were constructed from the high resolution crystal structure (1.91 Å) of a reference complex containing the inhibitor compound AHD1 (PDB code: 4FYR)39 using the graphical interface available in the molecular modeling program Insight-II40 and Discovery studio 2.541.
Molecular mechanics
Modeling of the AHD and PL ligand complexes was carried out by molecular mechanics using the CFF force field42 as described previously43.
Conformational research
The conformations of the free inhibitors were derived from their bound conformations in the PL complexes by gradual relaxation to the nearest local energy minimum, as previously described43.
Gibbs Free Energies Solvation
Ligand-receptor interactions take place in a solvent, which contributes to the binding process through hydrogen bonding and solvation phenomena. However, the electrostatic component of the Gibbs free energy (GFE) incorporating the effects of the ionic force through solving the nonlinear Poisson-Boltzmann equation44 was calculated by the Delphi module of Discovery Studio 2.541 as described previously43.
The calculation of binding affinity expressed as GFE complexation has been described in detail earlier43.
Interaction energy
The CFF force field was used to calculate the interaction energy (Eint) between the enzyme residues and the inhibitor, as previously reported43.
Generation of pharmacophores
Discovery Studio's 3D-QSAR (PH4) pharmacophore generation protocol41 via its Catalyst HypoGen algorithmic program45 was used to construct the APN inhibition PH4 as described previously43.
ADME properties
The pharmacokinetic profile of AHDs was calculated by the QikProp program46 as reported earlier43.
Virtual library generation
The generation of the virtual library was carried out as described in a previous study43.
ADME based library
The orientation of the virtual library was made using numerous selection criteria as described previously43.
Pharmacophore-based library search
The pharmacophore model (PH4) derived from the bound conformations of AHDs at the APN active site served as a library search tool, as previously described43.
Inhibitory power prediction
The conformer with the best mapping to the PH4 pharmacophore in each group of the targeted library subset was selected for in silico screening by the complexation QSAR model. The ∆∆Gcom calculation of each new selected analog was used to predict the APN inhibitory potency (Kipre) of the targeted AHD analog virtual library by inserting this parameter into the target-specific scoring function given in equation (1) parameterized using the AHD inhibitor training set complexation QSAR model39.
pKipre=−log10Kipre=a.∆∆Gcom + b (1)
RESULTS
Training and validation sets
Forty-six46 AHDs (Table 1) were selected from a series of compounds with experimentally determined properties and coming from the same laboratory39. Their experimental inhibitory activities (4.5 ≤ Kipre ≤ 4420 nM)39 cover a sufficiently wide range of concentrations to build a reliable QSAR model. The ratio between the sizes of the training and validation sets remains a critical point for correct classification but is limited by the number of sets of homologous compounds available in the literature47. In this study, a training set of 37 AHDs and a validation set of another 9 AHDs (Table 1) were created using the appropriate module of Discovery Studio 2.541.
One-descriptor QSAR model
Each of the 37 training sets (TS) and 9 validations sets (VS) APN: AHDx complexes (Table 1) was prepared by in situ modification of the crystal structure of the refined model (PDB entry code 4FYR) [39] of the APN: AHD1 complex as described in the Methods section. Additionally, the relative Gibbs free energy (GFE) of APN: AHDx ∆∆Gcom complex formation was calculated for each of the 46 optimized enzyme-inhibitor complexes. Table 1 lists the calculated values of ∆∆Gcom and its components as defined in equation (7), for the TS and VS of hydroxamic acid39. The QSAR model explained the variation of the experimental inhibitory potency of AHDs (pKiexp=– log10(Kiexp)) by correlating it with the GFE ∆∆Gcom calculated by linear regression (equation (1), Table 2), the validity of which by the statistical data of the regression is listed in Table 3, equation A and B. The correlation of ∆∆HMM and ∆∆Gcom explains approximately respectively 86 and 94 percent of the variation in pKiexp data and underlines the role of enthalpy contribution in ligand binding affinity.
The regression coefficient R2 of ∆∆Gcom attesting that structural information derived from 3D models of APN–AHDx complexes should conduct to reliable prediction of APN inhibitory potencies for novel AHD analogues (sharing the same binding mode) based on the QSAR B model, Table 3.
Binding mode of AHDs
The new series of Hydroxamic Acid used in this work has been synthesized39. Indeed, hydroxamic acids are used as metal ion chelators and the presence of the acid function in their molecular structure makes them particularly important for the inhibition of APN. Active site have been assessed from the X-rays crystal structure analysis of APN (PDB code 4FYR) in complex with one of the most active studied inhibitors in this work39.
Interaction Energy
The analysis of the interaction energy (IE) diagram per residue provides additional structural information to guide choice of the judicious R group to fill in the S1 and S1' pockets for AHD – APN binding affinity improvement. A comparative analysis of computed IE for the training set AHDs (Figure 4) divided into three classes (highest, moderate, and lowest activity) has been carried out to identify the residues for which the contribution to binding affinity could be increased.
However, the comparative study proves that we are the same level of IE contributions from active site residues for all three classes of inhibitors. Therefore, no suggestions of suitable substitutions able to improve the binding affinity as we previously reported for thymine-like inhibitors of APN. The statistical data confirmed validity of the correlation Equations (A) and (B) plotted on Figure 3. The ratio pKipre/pKiexp ≈1 (the Kipre values were estimated using correlation Equation B, Table 3) calculated for the validation set VAHD1-9 documents the substantial predictive power of the complexation QSAR model from Table 2. Thus, the regression Equation B (Table 3) and computed ∆∆Gcom GFEs can be used for prediction of inhibitory potencies Kipre against APN for novel AHD analogs, provided they share the same binding mode as the training set hydroxamic acid AHD1-37.
3D-QSAR Pharmacophore Model
Generation and validation of pharmacophore
The active conformation of 37 TS (AHD1-37) were used to generated APN inhibition 3D-QSAR pharmacophore and rated by 9 VS VAHD1-9 covering a large range of experimental activity (4.5 – 4420 nM) spanning almost three orders of magnitude. The three steps generation process: (i) the constructive, (ii) the subtractive, and (iii) the optimization step39 was described earlier43.
Hypotheses were scored according to errors in activity estimates from regression and complexity via a simulated annealing approach. The top scoring 10 unique pharmacophore hypotheses were kept, all displaying five-point features along with all relevant data listed in Table 4. They were selected based on significant statistical parameters, such as high correlation coefficient, low total cost, and low RMSD, Δ=3813.7= null cost (3863.3) – fixed cost (49.63); all meaning a high probability (>90%) that the model represents a true correlation43.
The evaluation of Hypo 1 is the mapping of the best active training set AHD1 (Figure 4 (D)) displaying the geometry of the Hypo1 pharmacophore of APN inhibition. The regression equation for pKiexp vs. pKipre estimated from Hypo1: pKiexp=1.0006× pKipre + 0.0028 almost equivalent to pKiexp=pKipre with a coefficient of 1 and in intercept of 0, plotted on Figure 4 (E) (n=37. R2=0.79. R2xv=0.78. F-test=130.03. σ=0.28, α > 95%) is also. Therefore, the PH4 is good potentially to choose the new AHD analogs.
Virtual Screening
In silico screening of a virtual (combinatorial) library can lead to hit identification as it was shown in our previous works on inhibitors design48-53.
Virtual Library
An initial Virtual library (VL) was generated by substitutions at positions for R1 and R2 (Table 5) on the scaffold. During the virtual library enumeration, the R-groups listed in Table 5 Were attached to in positions R1 and R2 of the AHD scaffold to form a combinatorial library of the size: R1×R2=252 252=56 644 analogs. All analogs are matching the substitution pattern of the best inhibitor AHD1. These AHDs analogs library was generated from fragments (chemicals) listed in databases of available chemicals54. To design a more target library of reduced size and increased content of drug-like molecules, we have introduced a set of filters and penalties such as the Lipinski rule-of-five (Mw>500 g/mol)55, which helped to select a smaller number of suitable AHDs that could be submitted to in silico screening.
In silico screening of library of AHDs
56 644 analogs of the library was further screening for molecular structures matching the 3D-PH4 pharmacophore model Hypo1 of APN inhibition. 95 AHDs mapped to at last 4 features of the pharmacophore according to the so-called similarity-property principle (SPP) according to which structurally similar compounds exhibit similar biological effects against the same target.
These best fitting analogs (PH4 hits) then underwent complexation QSAR model screening. The computed GFE of APN-AHDs complex formation, their components, and predicted inhibitory potency Kipre calculated from correlation Equation B (Table 3) is listed in Table 6).
Novel AHD analogs
The design of virtual library of novel analogs was guided by structural information retrieved from the AHDs active conformation and the pharmacophore model, used for the selection of appropriate substituents. The hydrophobic feature of PH4 at the position R1 shows clearly the type of group to be oriented towards the hydrophobic pocket S1. The analysis of frequency of occurrence of R-group during the selection of appropriate surrogates for two points of attachment: R1-group and R2- group shows that the frequency of occurrence of groups R1 and R2 among the best resulting from PH4 (Fig. 8) is as follows: for the large hydrophobic pocket S1’ filling R2-groups, 151: 4-(5-methyl-3-furyl)phenyl, 159: 4-(2-sulfanyl-cyclopenta-2,4-dien-1-yl)phenyl, 250: 2-sulfa-nylacetyl)amino, 115: 4-sulfanyl-5-(sulfanylmethyl)-pyrazol-1-yl, 57: 3-methanimidoylbenzenethiol, 242: 4-sulfamoylphenyl, 158: 4-[3,4-bis(sulfanyl) cyclo-penta-2,4-dien-1-yl]phe-nyl with occurrences of 16, 8, 5, 7, 7, 5 respectively are the most represented while 48: (3-sulfanylphenyl)for-mate with 2 occurrences, appears in the highest potency AHD analogs. In the smaller hydrophobic pocket S1, filling R1-groups 1: cyclopenta -2,4-diene-1-carbonyl, 17: amino (cyclopenta -2,4-dien-1-yl)methyl, 106: 5-fluoro-pyrazol-1-yl with occurrences of 5, 6, 5 is the most represented while 243: benzenesulfonyl and 126: pyrimidin-4-yl occurrences 1 and 2 appearing in the top 4 highest potency AHD analogs. The best analogs from these most commonly used substituents (R1-group: R2-group) are: 243- 242 (Kipre=0.05 nM); 126-242 (Kipre=0.07 nM); 126-158 (Kipre=0.09 nM); 17-48 (Kipre=0.13 nM). Branching larger aliphatic moieties in the R2 position for better filling the large S1’ pocket and conserving HB interactions and keeping almost the size in R1 position for the smaller hydrophobic pocket S1 of the AHD analogs contributed strongly to an overall improvement in the inhibitory activity against human M1 aminopeptidase (APN). This relates to the inhibitory potency of the best proposed new analogs.
Pharmacokinetic Profile of Novel AHD Analogs
The properties related to ADME such as Caco-2 cell permeability, blood-brain partition coefficient, octanol-water partitioning coefficient, aqueous solubility number of likely metabolic reactions, serum protein binding and another eighteen descriptors related to absorption, distribution, metabolism, and excretion (ADME) were calculated by the QikProp program46 for the new best AHD analogs (Table 7). The method of Jorgensen is used by this program56. Empirical data from more than 710 compounds including about 500 drugs and related hetero-cycles were used to produce regression equations correlating experimental and computed descriptors resulting in an accurate prediction of pharmacokinetic properties of molecules. Drug likeness (#stars) - the number of property descriptors that fall outside the range of optimal values determined for 95% of known drugs out of 24 selected descriptors computed by the QikProp, was used as an additional ADME-related compound selection criterion. The values for the predicted best active designed AHDs are compared to those computed for current anticancers targeting APN, displaying favorable pharmacokinetic profile with low number of stars indicating that the computed descriptors do not fall outside the range of 95% of known drugs (Table 7): (c) molar mass: 300 ≤ MW ≤ 500 g.mol-1; (d) total solvent-accessible molecular surface, (probe radius 1.4 Å): 300 ≤ Smol ≤ 1000 Å2; (e) hydrophobic portion of the solvent-accessible molecular surface, (probe radius 1.4 Å): 0 ≤ Smolhfo ≤ 750 Å2); (f) total volume of molecule enclosed by solvent-accessible molecular surface (probe radius 1.4 Å): 500 ≤ Vmol ≤ 2000 Å3; (g) number of non-trivial (not CX3), non-hindered (not alkene, amide, small ring) rotatable bonds: 0 ≤ #rotB ≤ 15; (h) estimated number of hydrogen bonds that would be donated by the solute to water molecules in an aqueous solution, values are averages taken over a number of configurations, so they can assume non-integer values: 0.0 ≤ HBdon ≤ 6.0; (i) estimated number of hydrogen bonds that would be accepted by the solute from water molecules in an aqueous solution, values are averages taken over a number of configurations, so they can assume non-integer values: 2.0 ≤ HBacc ≤ 20.0; (j) logarithm of partitioning coefficient between n-octanol and water phases: -2 ≤ LogPo/w ≤ 6.5; (k) logarithm of predicted aqueous solubility: logS, S in [mol.dm–3] is the concentration of the solute in a saturated solution that is in equilibrium with the crystalline solid: -6.0 ≤ LogSwat ≤ 0.5; (l) logarithm of predicted binding constant to human serum albumin: -1.5 ≤ LogKhsa ≤ 1.5; (m) logarithm of predicted brain/blood partition coefficient: -3.0 ≤ LogB/B ≤ 1.2; (n) predicted apparent Caco-2 cell membrane permeability in Boehringer-Ingelheim scale in [nm s-1]: BIPcaco < 25 poor, BIPcaco > 500 nm.s-1 great; (o) number of likely metabolic reactions: 1 ≤ #metab ≤ 8; (p) predicted inhibition constants Kipre. Kipre in pM was predicted from computed ∆∆Gcom using the regression Equation B shown in Table 3; (q) HOA: human oral absorption: 1=low, 2 =medium, 3=high; (r) % HOA: percentage of human oral absorption in gastrointestinal tract: ≥ 80%=high.
DISCUSSION
The most comprehensive metrics of APN inhibition by hydroxamic acid containing AHDs reported by J. Lee et al.39. Intermolecular interactions of AHD1 and hAPN including hydrophobic stacking interactions and hydrogen bonds were the key determinants for better affinity with the target. The exploration of the chemical AHD subspace implemented in a diverse virtual library with AHDs active conformation yielded the best R1 and R2 substituent’s to be accommodated by the hydrophobic pockets or rooted in other ways such as hydrogen bonds and van der Waals contacts. The strategy was executed over three orders of magnitude of experimental Ki, i.e. three pKi units taking benefit from the reported SAR continuity39 making feasible activity prediction according similarity-property principle (SPP).
The compound 6f, N-(2-(Hydroxyamino)-2-oxo-1-(3′-fluoro-[1,1′-biphenyl]-4-yl)ethyl)-4-(methyl-sulfon-amido) benzamide has been designed by J. Lee et al., with the purpose to improve both potency and solubility through removal of two fluorine atoms to keep only one compared to AHD1 (Ki=4.5±0.8 nM), they reached a potency Ki=0.66 ± 0.06 nM57. Used AHD analogs potency prediction model computed ∆∆Gcom=- 2 kcal/mol and a potency Ki=2.1 nM using correlation Equation B, Tables 3 and 6, presenting 6f as twice more potent than AHD1 and keeping in this way the same trend as experimental values according to which, 6f is 6-fold more potent than AHD1 regardless experimental uncertainties. The computed solubility of some AHD analogs (Table 7) is of the same order as of 6f.
The predicted most potent analogs 243-242 (50 pM) with benzenesulfonyl (243) in R1 and 4-sulfamoylphenyl (242) in R2, 126-242 (70 pM) bearing pyrimidin-4-yl in R1, 126-158 (90 pM) with 4-[3,4-bis(sulfanyl)cyclopenta-2,4-dien-1-yl]phenyl (158) in R2 keep the filling of S1 bringing better interactions and fill better the large S1’ hydrophobic pocket resulting in better affinity as displayed in Figure 10 comparing the interaction energy breakdown to APN active site residues of the best active TS AHD1 and novel analogs. This substantial stabilization will undergo medicinal chemistry verification through synthesis and biological evaluation.
Limitations of the study
The main limitation of this MM – PB study is the lack of experimental verification of the predicted novel analogs potency. Nevertheless the novel AHD analogs – APN complexes’ stabilization is cross checkable through Molecular Dynamics runs in order to confirm the active site residues’ side chains stabilized conformation and by the way that of the novel most potent AHD analogs as presented (see Figure 8 for example). Usually this last check is a relevant step before synthesis and biological evaluation for those who cannot afford assays step easily. Unfortunately these time-consuming MD runs represent a tremendous effort we’re preparing to address in due course.
CONCLUSION
SAR structural investigation of hydroxamic acid derivatives as a novel human M1 aminopeptidase (APN) cancer inhibitor from the crystal structure of APN: AHD complex guided us while preparing a QSAR model for the reliable complexation of APN activation that correlates with the calculated relative Gibbs free energies to form a complex with observed APN activation potencies. In addition we have derived a 3D-QSAR PH4 model from AHD active conformation using a training set of 37 and validation set of 9 AHDs with known activation activities. Careful analysis of interactions between the APN’s active site residues and APNs directed us in the design of an initial diversity virtual combinatorial library of new AHD analogs with multiple substitutions of hydrophobic groups in R1 and R2. A library screened by matching of the analogs to the PH4 pharmacophore permitted selection of a library subset of AHDs. This subset of 95 best virtual hits was submitted to computation of predicted activation potencies by the complexation QSAR model. The hit analogs reached predicted activities in the picomolar concentration range. The hit designed AHD analogs 243-242 (50 pM), 126-242 (70 pM), 126-158 (90 pM) and 17-48 (130pM) are recommended for synthesis and subsequent activity evaluation in APN activation assays and may lead to a discovery of novel hydroxamic potent partial APN agonists.
ACKNOWLEDGEMENTS
Authors are thankful for Fundamental and Applied Physics Laboratory, University Nangui Abrogoua, Ivory Coast to provide necessary facilities for this work.
AUTHOR’S CONTRIBUTIONS
This work was carried out in collaboration among all authors. All authors read and approved the final manuscript.
COMPETING INTERESTS
The authors declare there is no conflict of interest in relation with the work presented herein.
REFERENCES