A Single Point Mutation Blocks the Entrance of Ligands to the Cannabinoid CB2 Receptor via the Lipid Bilayer
Laboratory of Computational Medicine, Biostatistics Unit, Faculty of Medicine, Universitat Autònoma Barcelona, 08193 Bellaterra, Barcelona, Spain
Department of Biochemistry and Physiology, Faculty of Pharmacy and Food Sciences, Universitat de Barcelona, 08028 Barcelona, Spain
Centro de Investigación en Red, Enfermedades Neurodegenerativas (CIBERNED), Instituto de Salud Carlos III, 28031 Madrid, Spain
Department of Biochemistry and Molecular Biomedicine, Faculty of Biology, Universitat de Barcelona, 08028 Barcelona, Spain
*E-mail: leonardo.pardo@uab.cat.Abstract
Molecular dynamic (MD) simulations have become a common tool to study the pathway of ligand entry to the orthosteric binding site of G protein-coupled receptors. Here, we have combined MD simulations and site-directed mutagenesis to study the binding process of the potent JWH-133 agonist to the cannabinoid CB2 receptor (CB2R). In CB2R, the N-terminus and extracellular loop 2 fold over the ligand binding pocket, blocking access to the binding cavity from the extracellular environment. We, thus, hypothesized that the binding pathway is a multistage process consisting of the hydrophobic ligand diffusing in the lipid bilayer to contact a lipid-facing vestibule, from which the ligand enters an allosteric site inside the transmembrane bundle through a tunnel formed between TMs 1 and 7 and finally moving from the allosteric to the orthosteric binding cavity. This pathway was experimentally validated by the Ala2827.36Phe mutation that blocks the entrance of the ligand, as JWH-133 was not able to decrease the forskolin-induced cAMP levels in cells expressing the mutant receptor. This proposed ligand entry pathway defines transient binding sites that are potential cavities for the design of synthetic modulators.
1Introduction
The approximately 800 (∼450 for sensory and ∼350 for nonsensory functions) members of the G protein-coupled receptor (GPCR) family are mainly formed by a conserved architecture of seven transmembrane domains (TMs).1 The intracellular domains and/or the intracellular part of the TM segments bind a small repertoire of signaling proteins (16 G proteins, four arrestins, or seven G protein-coupled receptor kinases),2−4 whereas the extracellular domains (ECDs) and/or the extracellular part of the TM segments recognize exogenous signals (e.g., 1 trillion olfactory5 or ∼1 000 bitter molecules6), or endogenous signals that range in size from ions (e.g., Ca2+) to small molecules (e.g., glutamate, neurotransmitters) to peptides (e.g., opioids, endothelin, glucagon) or large proteins (e.g., chemokines, glycoprotein hormones), or ∼200 000 synthetic ligands,7 or ∼475 drugs.8 Thus, the mechanism of binding of this set of highly diverse extracellular ligands varies considerably depending on the ligand and receptor,9 as subsequently confirmed by the analysis of known inactive and active structures of GPCRs.10
The high druggability of membrane embedded GPCRs11 is due to the fact that their malfunction, with 435 disease-associated mutations,12 commonly translates into pathological outcomes13 and because the orthosteric binding site, a conserved pocket within the 7TM domain that optimally accommodates the electrostatic and steric properties of the ligand, is normally accessible from the extracellular space. However, the depth of ligand penetration into the TM bundle is considerable,14 so ligands must also engage stable interactions at the entrance of the binding site. Ligand binding pathway simulations on β-adrenergic receptors have identified these additional cavities that are transiently occupied before arrival to the orthosteric binding site.15,16 They have been named extracellular vestibules15 or entrances,17 or secondary16,18 or metastable19 binding sites, or exosites.20 Cavities like those described for β-adrenergic receptors have also been found for the M3 muscarinic receptor,21 the adenosine A2A receptor,22 the histamine H4 receptor,23 the dopamine D2 and D3 receptors,24 and olfactory receptors,25 among others. Remarkably, while the orthosteric binding site is conserved among family members, these sites at the ECDs are highly divergent and have become potential binding sites for synthetic modulators.26 Thus, they have been used in the design of allosteric modulators27 or bitopic ligands,20,28 as well as proposed to bind short bivalent ligands.29,30
Class A GPCRs activated by hormone-like signaling molecules derived from lipid species with long hydrophobic moieties31 possess a distinctive structural signature at the ECD, in comparison with other class A GPCRs that are activated by polar ligands.32 In the crystal structures of S1P1,33 LPA1,34 FFAR1,35 CRTH2,36 CB1,37,38 and CB239 receptors, the extracellular N-terminus and extracellular loop 2 fold over the ligand binding pocket. Thus, the entrance of ligands to the orthosteric binding site through the membrane bilayer has been proposed for S1P1,33,40 CRTH2,36 and CB137,41 receptors via TMs 1 and 7, for the FFAR1 via TMs 3 and 4,35 for the MT1 melatonin receptor via TMs 4 and 5,42 for opsin via TMs 5 and 6,43,44 and for CB245 and PAR146 receptors via TMs 6 and 7. In contrast, some authors propose that hydrophobic ligands of LPA134 and CB1 and CB247 receptors may reach the orthosteric site from the extracellular environment.
In this manuscript, we combine molecular dynamic (MD) simulations and site-directed mutagenesis to study the binding process of a potent and selective CB2R agonist, JWH-133. The proposed binding pathway is a multistage process consisting of the ligand diffusing in the lipid bilayer to contact a lipid-facing vestibule, from which the ligand enters an allosteric site inside the TM bundle through a tunnel formed between TMs 1 and 7, and finally moving from the allosteric to the orthosteric binding cavity.
2Materials and Methods
2.1Initial CB2R Models
The inactive structure of CB2R was used (PDB id 5ZTY).39 Protonation states were assigned with PROPKA.48 The system was oriented by the Orientations of Proteins in Membranes (OPM) database49 and embedded in four different lipid bilayer boxes (9 × 9 × 9.5 nm), constructed using PACKMOL-memgen,50 containing 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), cholesterol (CHL; 10:6 POPC/CHL ratio), water molecules (TIP3P), and monatomic Na+ and Cl– ions (0.15 M). JWH-133 (JWH) was inserted into the membrane by replacing one in six CHL molecules from both leaflets, leading to a 10:5:1 POPC/CHL/JWH ratio. This substitution of CHL by JWH was randomly performed five times for each of the four constructed membrane systems (see Figure S1).
2.2Molecular Dynamics Simulations
These 20 combinations of CB2R models were simulated with GROMACS2016.4.51 The amber14sb-ildn force field was used for the protein, solvent. and ions;52 a GROMACS adaptation of lipid14 for lipids;53,54 and the general Amber force field (GAFF2) with HF/6-31G*-derived RESP atomic charges for JWH-133.55 Molecular systems were subjected to 10 000 steps of energy minimization, followed by 20 ns of gradual relaxation of positional restraints (corresponding to 100, 50, 25, and 10 kJ mol–1 nm–2) at protein backbone coordinates before the production phase to hydrate the cavities and allow lipids to pack around the protein. After equilibration, five replicas of 1 μs unrestrained MD trajectory, of each of the 20 combinations, was generated at a constant temperature of 300 K using separate v-rescale thermostats for the receptor, ligand, lipids, and solvent molecules. A time step of 2.0 fs was used for the integration of equations of motions. Bonds involving hydrogen atoms were kept frozen using the LINCS algorithm. Lennard-Jones interactions were computed using a cutoff of 1.1 nm, and the electrostatic interactions were treated using PME with the same real-space cutoff under periodic boundary conditions (see Figure S1). The analysis of the trajectories was performed using MDAnalysis;56 the stability of the membrane was monitored (see Figure S2) using FATSLiM.57 Visualization and image rendering were performed with PyMOL58 and VMD,59 and the tunnels of CB2R, from the orthosteric binding cavity, were calculated with the CAVER program.60
2.3Metadynamics Simulations
Metadynamics simulations (see Figure S1) were performed using GROMACS2016.4 patched with PLUMED version 2.561.62 The starting point was taken from one of the 100 unbiased MD trajectories in which JWH-133 spontaneously bound the lipid-facing part of TMs 1 and 7 (orange in Figure 1b). The end point is the docking pose of JWH-133 to the orthosteric binding site of CB2R previously reported,63 using the structurally similar AM12033 ligand (PDB id 6KPC)47 as a template. The distance between the center of mass at the starting and end points of JWH-133 was used as a collective variable (CV). This CV was selected for simplicity and because it can efficiently represent the dimensional space of the free-energy of binding. In metadynamics simulations, a history-dependent potential is applied along the CV, built as a sum of Gaussian kernels, with a specific width of 0.1, height of 0.48 kJ/mol, and pace of 5000 steps (i.e., 0.1 ns). A light RMSD constraint (with an energy constant KAPPA of 200 IU) was applied to ensure a correct end conformation of JWH-133 at the orthosteric site. We used well-tempered metadynamics, where the height of the applied Gaussian potentials is rescaled so that the system does not explore free-energy regions beyond the expected values. Thus, a bias factor of 20, accounting for the difference between the system temperature and the temperature of the CV, was set up. Convergence of the well-tempered metadynamics simulations was determined when JWH-133 had explored all of the CV space and the free-energy profile was constant at 10 ns intervals. Free-energy profiles were calculated using PLUMED sum_hills function.
2.4CB2R Mutants
Mutations were produced using the QuikChange Site-Directed Mutagenesis Kit. The cDNA for hCB2R, cloned into pcDNA3.1, was amplified using sense and antisense primers harboring the triplets for the desired point mutation (Pfu turbo polymerase was used). The forward and reverse primers are ccccagaagacagctTTTgctgtgttgtgcactc and gagtgcacaacacagcAAAagctgtcttctgggg for Val361.35Phe, ccccagaagacagctATGgctgtgttgtgcactc and gagtgcacaacacagcCATagctgtcttctgggg for Val361.35 Met, gtcaagaaggcctttTTCttctgctccatgctg and cagcatggagcagaaGAAaaaggccttcttgac for Ala2827.36Phe, and gtcaagaaggcctttATGtctgctccatgctg and cagcatggagcagaaCATaaaggccttcttgac for Ala2827.36Phe, respectively, in which mutated nucleotides are written in upper case characters. The nonmutated DNA template was digested for 1 h with DpnI. PCR products were used to transform XL1-blue supercompetent cells. Finally, positive colonies were tested by sequencing to select those expressing the correct DNA sequence.
2.5cAMP Determination Assays
Determination of cAMP levels in HEK-293T cells transiently expressing CB2R (1 μg of cDNA) or the mutant version of the receptor was performed using the Lance-Ultra cAMP kit (PerkinElmer). Two hours before initiating the experiment, the medium was substituted by a serum-free medium. Then, transfected cells were dispensed in white 384-well microplates at a density of 3000 cells per well and incubated for 15 min at room temperature with compounds, followed by 15 min of incubation with forskolin and 1 h more with homogeneous time-resolved fluorescence (HTRF) assay reagents. Fluorescence at 665 nm was analyzed on a PHERAstar Flagship microplate reader equipped with an HTRF optical module (BMG Labtech). Data analysis was made based on the fluorescence ratio emitted by the labeled cAMP probe (665 nm) over the light emitted by the europium cryptate-labeled anti-cAMP antibody (620 nm). A standard curve was used to calculate cAMP concentration. Forskolin-stimulated cAMP levels were normalized to 100%.
3Results
3.1Computational Model to Study the Pathway of Ligand Entry to CB2R
The selected ligand is JWH-133, a potent selective CB2R agonist.64 We have chosen JWH-133 over other agonists for its rigidity and for having only one heteroatom in its structure. JWH-133 (6.4) and cholesterol (8.7) have similar computed logP values, which makes JWH-133 a potential through-the-membrane binder. The model to study the entry of JWH-133 into the orthosteric binding site consisted of inactive CB2R embedded in a lipid bilayer box (see Materials and Methods). Four different membranes were constructed, and for each of them, JWH-133 was positioned at five different positions. These 20 combinations were each subjected to five rounds of unbiased 1 μs MD simulation with an aggregate sampling of 100 μs (see Methods and Figure S1). Data was collected every 10 ns.
3.2Tunnels in CB2R
To delineate the pathways from the orthosteric binding cavity of CB2R, buried in the TM bundle, to the surrounding solvent, either the aqueous extracellular environment or the lipid bilayer, we used the ensemble of structures collected during the MD simulations of CB2R. Figure 1a displays three different clusters of tunnels calculated with the skeleton search algorithms implemented in the CAVER program.60 These are toward the extracellular domain and toward the lipid bilayer via either TMs 2 and 3 or TMs 1 and 7. The tunnel route preference was evaluated using throughput, which is the probability of ligand transportation throughout the pathway (the higher the value, the greater the importance of the pathway). These values of throughput (Figure 1a), which consider the conformational flexibility of CB2R, suggest that the transit of molecules through the lipidic phase via TMs 1 and 7 is favored compared to TMs 2 and 3 or to the extracellular route.
3.3Ligand Diffusion to CB2R and Tunnel Opening
None of the JWH-133 agonists, embedded at different positions in the membrane, spontaneously bound the CB2R orthosteric binding site during the aggregated 100 μs of MD simulation, but in three simulations the ligand remained bound to the receptor surface, mainly to amino acids in TMs 5 and 6 (one simulation) or TMs 1 and 7 (two simulations). Figure 1b shows these stable surface-bound ligand conformations with root-mean-square deviation (rmsd) values < 20 Å relative to the reference docked binding mode of JWH-133 (see Materials and Methods), while they have displaced >20 Å from their initial position. These results, together with the throughput values of tunnel route preference (see above), reinforce TMs 1 and 7 as the most feasible pathway of ligand entry to the orthosteric site of CB2R. Figure 1c shows the key Val361.35, Leu391.38, Phe2837.37, and Met2867.40 residues (Ballesteros and Weinstein numbering scheme65 as superscript) located at the entrance of TMs 1 and 7. Of note is the conformation of Phe2837.37, which in the gauche+ conformation, like that observed in the crystal structures of CB2R,39,47,66 closes the tunnel, whereas in the trans conformation, it opens it. Moreover, it has been shown that the N-terminus of the S1P1 receptor packs against ECL-2 in the active conformation, leading to an opening of the ligand access vestibule between TMs 1 and 7.67 Thus, we have quantified in Figure S3 the number of snapshots, during the 100 μs of MD simulation, in which Phe2837.37 adopts the trans conformation or/and the distance between the top of TMs 1 (Cα atom of Thr341.33) and 7 (Cα atom of Val2777.31) increases from the initial value of 14.5 Å to values larger than 15.5 Å. GPCRs are dynamic proteins that permit rapid small-scale structural fluctuations,68 and accordingly, both events can be freely observed, simultaneously or not, during the calculated trajectories (Figure S3). Importantly, the trajectory in which the ligand (orange in Figure 1b) spontaneously binds the lipid-facing part of TMs 1 and 7 (surface-bound, Figure S3) triggers or stabilizes the trans conformation of Phe2837.37 (Figure 1d), and the distance between TMs 1 and 7 increases to 16.5 Å (Figure 1f) most of the simulation time. In contrast, in other trajectories in which the ligand did not bind the receptor surface and these events did not spontaneously occur (surface-unbound, Figure S3), Phe2837.37 remains in the gauche+ conformation (Figure 1d) and the distance between TMs 1 and 7 decreases to values lower than the initial 14.5 Å (Figure 1f). Figure 1e shows two selected snapshots of these surface-unbound and surface-bound trajectories (see Figure S3) in which the tunnel between TMs 1 and 7 is closed and open, respectively.
3.4Ligand Access to the Orthosteric Binding Site
The above results, and the experimental evidence (see Introduction), led us to use the ligand bound to TMs 1 and 7 of CB2R to progressively sample the binding event. Because entry to the orthosteric binding site can have the largest energetic barrier of the process,15,40 in some cases unsurmountable by unbiased MD simulations, we used well-tempered metadynamics.69 In this technique, a biasing potential is applied to permit the system to explore energetically inaccessible regions for unbiased MD simulations on reasonable time scales. We used as a collective variable of the distance between the center of mass of the ligand in the initial conformation and the reference docking pose at the orthosteric binding site (see Materials and Methods). In the free-energy profile of five independent simulations (Figure 2a and b), three minima can be noticed in most trajectories. These steps in the entry pathway were designated as lipid-facing, bundle-facing, and orthosteric. In the first minimum, the ligand is bound to the lipid-facing part of the CB2R tunnel mainly formed by amino acids (Val361.35, Cys401.39, Leu431.42, Ala2827.36, Phe2837.37, Met2867.40, Leu2897.43, Ile2907.44, and Met2937.47) in TMs 1 and 7 (Figure 2c). In the second minimum, the ligand is bound to the bundle-facing part of the CB2R tunnel mainly formed by amino acids (Val361.35, Cys401.39, Phe872.57, Phe912.61, His952.65, Ala2827.36, and Met2867.40) in TMs 1, 2, and 7 (Figure 2d). In the final third minimum of the binding process, JWH-133 reaches the orthosteric binding site. This computationally derived binding pose, using metadynamics, reproduced the main contacts between JWH-133 and the receptor (Phe872.57, Phe912.61, Phe942.64, His952.65, Thr1143.33, Phe1173.36, Phe183ECL2, Trp1945.43, Phe2817.35, and Ser2857.39) previously proposed63 (Figure 2e).
3.5Experimental Validation of the Pathway of Ligand Entry
To experimentally validate the proposed mechanism of ligand entry by site-directed mutagenesis, we explored the amino acids involved in the transition from the lipid- to the bundle-facing part of the CB2R tunnel (Figure 2c,d). Val361.35 and Ala2827.36 are located midway between these two minima, delimiting the tunnel between TMs 1 and 7. Notably, three aromatic, Phe912.61, Phe942.64, and His952.65, side chains in TM 2 are in their environment (Figure 3a). Thus, we obtained the Val361.35Phe and Ala2827.36Phe mutant versions of the CB2R with the main idea that the aromatic side chains would form a stable aromatic cluster with Phe912.61. Phe942.64, and/or His952.65, blocking the access of the ligand to the orthosteric binding site. Dose–response experiments in HEK-293T cells expressing the CB2R stimulated with forskolin and treated with JWH-133 showed reduced cAMP production (Figure 3b and Table 1), as expected for Gi-coupled receptors. Clearly, both Val361.35Phe and Ala2827.36Phe mutations impair the effect of JWH-133, as the decrease of forskolin-induced cAMP is of less magnitude in the mutant receptors (the effect is more noticeable in Ala2827.36Phe than in Val361.35Phe).
The importance of these stable aromatic interactions was tested by comparing the results with those obtained using Val361.35 Met and Ala2827.36 Met mutants (Figure 3b and Table 1). The highly polarizable sulfur atom of the Met side chain can form with other aromatic side chains stronger interactions than Phe.70,71 Interestingly, these new mutations do not affect the decrease in forskolin-induced cAMP in a significant manner, thus suggesting they do not block the entrance of the JWH-133 ligand probably due to the higher flexibility of Met relative to Phe. Our experimental results point out that the highly stable aromatic cluster between Phe912.61, Phe942.64, His952.65, and Phe2827.36 (in the Ala2827.36Phe mutation) blocks the entrance of JWH-133 to the orthosteric binding site.
3.6Ligand Access to the Binding Site in the Val361.35Phe and Ala2827.36Phe Mutant Receptors
To study the influence of this aromatic cluster, attained in the Val361.35Phe and Ala2827.36Phe mutants, during the process of ligand entry, we performed well-tempered metadynamics using the same protocol as for nonmutated CB2R. As expected, the Val361.35Phe mutation hampers the ability of the ligand to access the orthosteric site in one of the replicates and the possibility to explore the minimum corresponding to the bundle-facing part of the CB2R tunnel in three other replicates (Figure 3c). On the other hand, the Ala2827.36Phe mutation blocks the entrance of the ligand in four of the five replicates performed (Figure 3d) in agreement with our hypothesis and experimental results.
4Discussion and Conclusions
Drug-target residence time is a key element as a biological effect requires binding to the target receptor.72,73 Residence time is related to the ligand–receptor binding kinetics and, thus, to the pathway used in ligand association and dissociation.74 MD simulations have become a convenient tool to estimate residence times75,76 and to analyze such pathways.15,16,21−25,40,77 These simulations show that ligand binding is a multistage process in which the ligand forms key interactions with amino acids at the entrance of the binding site. Most of these simulations are for class A GPCRs that bind polar ligands directly from the extracellular environment.
Here, we have studied the binding pathway of the potent JWH-133 agonist to CB2R by MD simulations. In CB2R, the N-terminus and ECL-2 that lacks the disulfide bridge to TM 3, characteristic of the GPCR family, fold over the ligand binding pocket39 blocking access to the orthosteric binding cavity from the extracellular environment. We, thus, propose that the hydrophobic JWH-133 (calculated logP of 6.4) diffuses through the bilayer leaflet to contact a lipid-facing cavity of CB2R mainly formed by amino acids in TMs 1 and 7 (Figure 4). In this point of the pathway, the ligand forms hydrophobic/aromatic contacts with Leu391.38, Leu421.41, Leu431.44, Phe2837.37, Met2867.40, and Ile2907.44 (Figure 2c). Importantly, the binding of the ligand to this lipid-facing cavity triggers and/or stabilizes the trans conformation of Phe2837.37 (Figure 1d) and increases the distance between TMs 1 and 7 (Figure 1f), opening the tunnel between TMs 1 and 7 (Figure 1e).
Subsequently, JWH-133 moves to a second energetic minimum in the pathway, located in the bundle-facing part of the CB2R tunnel (Figure 4), to mainly interact with Val361.35, Cys401.39, Phe872.57, Phe912.61, His952.65, Ala2827.36, and Met2867.40 (Figure 2d). We experimentally validated this entrance tunnel by mutating Ala2827.36 that is not located at the orthosteric binding site but midway between the lipid- and the bundle-facing part of the CB2R tunnel (Figure 3b). It is common in the GPCR field to perform mutagenesis assays to study the influence in affinity or efficacy of amino acids located at the orthosteric binding pocket. Here, we aimed at decreasing the residence time of the ligand by blocking its entrance and, thus, blocking the signaling response (decrease of forskolin-induced cAMP). The Ala2827.36Phe mutation, but not Ala2827.36 Met, allows a blockade of the effect of JWH-133 because the decrease of forskolin-induced cAMP is of less magnitude in the mutant than in nonmutated receptor (Figure 3b). We propose that a highly stable aromatic cluster between Phe912.61, Phe942.64, His952.65, and Phe2827.36 (in the Ala2827.36Phe mutation), which cannot be accomplished with Met2827.36, blocks the entrance of JWH-133 to the orthosteric binding site. Figure S4 shows a structure-based sequence alignment of the amino acids forming this pathway among the members of class A GPCRs that bind hormone-like signaling molecules derived from lipid species. Finally, JWH-133 moves to the most energetically stable position in the pathway, at the orthosteric binding site (Figure 4).
The pathway of JWH-133 entry to CB2R defines two transient binding sites, in addition to the final location at the orthosteric cavity (Figure 4). These transient sites have been proposed to be potential binding sites for synthetic modulators (allosteric modulators or bitopic ligands), becoming important sites for drug discovery due to the high conservation of the traditional orthosteric binding site.18,19,79 To analyze this suggestion, we have superimposed our previously reported binding modes of a negative allosteric modulator63 and a bitopic ligand78 of CB2R. Clearly, the bitopic ligand binds at the orthosteric site and the identified transient lipid-facing site. Thus, the lipid-facing site in CB2R is analogous to the vestibule or secondary or metastable binding site defined for other GPCRs that bind polar ligands (see Introduction), which is key for ligand desolvation15,16,40 and for selectivity,20,28,77 but it might not allosterically modulate receptor activity. In contrast, the negative allosteric modulator cannabidiol binds at the identified bundle-facing site, from which it elicits the allosteric modulation.63 Thus, the bundle-facing site is an allosteric site in the CB2R, located near the orthosteric site, as has been suggested for muscarinic receptors.79 In conclusion, our findings have shown that ligand pathway simulations might be a useful tool to identify and delineate cavities for the design of synthetic modulators of GPCRs.
5Data and Software Availability
PACKMOL-Memgen, distributed with AmberTools, is free of charge. MDAnalysis, GROMACS, and PLUMED are open source. FATSLiM and CAVER are licensed under the GNU General Public License. VMD is available to noncommercial users under a distribution-specific license, and PyMOL is commercial software with a paid license.
Supporting Information Available
The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jcim.2c00865.
- Workflow of the computational process, heatmap representation of the cavity opening events, sequence alignment of class A GPCRs that bind hydrophobic-like signaling molecules, and analysis of membrane stability (PDF)
Supplementary Material
Notes
The authors declare no competing financial interest.
Acknowledgments
We acknowledge the financial support from the Spanish Ministry of Economy and Innovation with FEDER funds (RTI2018-094204-B-I00, PID2019-109240RB-I00, PID2020-113430RB-I00). C.L.T. is recipient of a FPI fellowship (BES-2017-081872), and M.G.-A. and A.Q.G. acknowledge the Universitat Autònoma de Barcelona for their predoctoral grant.