- Split View
-
Views
-
Cite
Cite
Louis Gauthier, Rémicia Di Franco, Adrian W R Serohijos, SodaPop: a forward simulation suite for the evolutionary dynamics of asexual populations on protein fitness landscapes, Bioinformatics, Volume 35, Issue 20, October 2019, Pages 4053–4062, https://doi.org/10.1093/bioinformatics/btz175
- Share Icon Share
Abstract
Protein evolution is determined by forces at multiple levels of biological organization. Random mutations have an immediate effect on the biophysical properties, structure and function of proteins. These same mutations also affect the fitness of the organism. However, the evolutionary fate of mutations, whether they succeed to fixation or are purged, also depends on population size and dynamics. There is an emerging interest, both theoretically and experimentally, to integrate these two factors in protein evolution. Although there are several tools available for simulating protein evolution, most of them focus on either the biophysical or the population-level determinants, but not both. Hence, there is a need for a publicly available computational tool to explore both the effects of protein biophysics and population dynamics on protein evolution.
To address this need, we developed SodaPop, a computational suite to simulate protein evolution in the context of the population dynamics of asexual populations. SodaPop accepts as input several fitness landscapes based on protein biochemistry or other user-defined fitness functions. The user can also provide as input experimental fitness landscapes derived from deep mutational scanning approaches or theoretical landscapes derived from physical force field estimates. Here, we demonstrate the broad utility of SodaPop with different applications describing the interplay of selection for protein properties and population dynamics. SodaPop is designed such that population geneticists can explore the influence of protein biochemistry on patterns of genetic variation, and that biochemists and biophysicists can explore the role of population size and demography on protein evolution.
Source code and binaries are freely available at https://github.com/louisgt/SodaPop under the GNU GPLv3 license. The software is implemented in C++ and supported on Linux, Mac OS/X and Windows.
Supplementary data are available at Bioinformatics online.
1 Introduction
Protein coding sequence evolution broadly depends on two questions. First, how random mutations change protein structure, function, and consequently organismal fitness. Second, which of these random mutations eventually survive or are purged in evolution. The first question relates to the distribution of fitness effects of random mutations (DFE), which can be determined from the knowledge of the fitness landscape or genotype-phenotype relationship. The second question relates to the role of population size, dynamics and structure in determining a mutation’s fixation probability. Unfortunately, these two questions seldom intersect—the first is traditionally asked by molecular biochemists and biophysicists, and the second is asked by evolutionary biologists and population geneticists. Despite the efforts to combine these two causalities in molecular evolution (Bershtein et al., 2017; DePristo et al., 2005; Echave and Wilke, 2017; Goldstein, 2011; Harms and Thornton, 2013; Liberles et al., 2012; Silander et al., 2007), there remains a divide between these disciplines, both in concepts and methods. Indeed, to date, there is no publicly available computational tool that integrates molecular biophysics and population genetics.
Currently available methods to perform forward evolutionary simulations differ in scale, scope and flexibility, but they are generally intended to investigate potential scenarios for whole genome evolution based on observed genetic variation among natural populations (Carvajal-Rodriguez, 2008; Hoban et al., 2012). Several forward-based tools including forqs (Kessner and Novembre, 2014), SLiM (Messer, 2013), ForwSim (Padhukasahasram et al., 2008), GENOMEPOP (Carvajal-Rodriguez, 2008), FFPopSim (Zanini and Neher, 2012), QuantiNemo (Neuenschwander et al., 2008), GeneEvolve (Tahmasbi and Keller, 2017), simuPOP (Peng and Kimmel, 2005), SFS_CODE (Hernandez, 2008) and fwdpp (Thornton, 2014) implement genetic features such as chromosome types, linkage and recombination. However, those tools concern themselves with patterns of polymorphism and structural variation across chromosomes, rather than explicit coding sequence evolution. Thus, it is challenging to model the evolution of protein sequences. Other programs such as OncoSimulR (Diaz-Uriarte, 2017) model the evolution of large asexual populations with user-defined fitness landscapes but enforce strictly bi-allelic loci on limited sites and do not model DNA sequences explicitly. Altogether, these tools account for neither the biochemical nor biophysical features of specific gene products.
Another class of models in protein evolution are the methods in molecular phylogenetics. Embedded in these phylogenetic approaches is a quantitative model of protein evolution that describes the amino acid substitution rates (Rodrigue et al., 2010). Substitution matrices contain information on the rate at which mutations can arise and the rate at which they can fix based on their estimated selective advantage (Yang and Nielsen, 2002). These transition matrices may also contain implicit information on the biophysics of proteins—for instance, the transition probabilities between amino acids of similar chemical properties are higher than those amino acids of different types. In some cases, these matrices can also include information on the tertiary structure of proteins (Halpern and Bruno, 1998; Lartillot and Philippe, 2004; Scherrer et al., 2012). However, none of these models explicitly account for the contribution of population dynamics and structure in shaping sequence evolution.
Finally, there are the physics-based models of protein evolution. To investigate the role of biophysical properties and structure to protein evolution (Bloom et al., 2007; Shakhnovich, 2006; Taverna and Goldstein, 2000), these models often rely on simplified representations of proteins that fold their sequences on a lattice to calculate biophysical properties. As such, they can estimate folding stability and protein-protein interactions. Nonetheless, most of these models for protein evolution are agnostic to the role of population size and dynamics. Although there are studies that investigated the interplay between population structure and protein biophysics (Rotem et al., 2018; Serohijos et al., 2013; Wylie and Shakhnovich, 2011), the computational tools for performing forward evolutionary simulations that account for both factors are not yet available for the community.
More broadly, this synthesis of population genetics and protein biophysics is important in the evolution of microbes and pathogens, such as the acquisition of antibiotic resistance and viral evolution. Several works have used this integrative approach to explore the interplay between different scales in evolutionary biology. For example, Rotem et al. showed that the evolution and dynamics of biophysical traits of an RNA virus subjected to a neutralizing antibody are strongly dependent on population size (Rotem et al., 2018). Salverda et al. showed how the dynamics of TEM-1 β-lactamase adaptation is determined by the topography of the fitness landscape and by mutational supply (Salverda et al., 2017). Finally, Heckmann et al. used a population-genetic model to show that enzyme kinetic parameter evolution in Escherichia coli is constrained by strong epistatic interactions (Heckmann et al., 2018).
Here we introduce SodaPop, an efficient forward simulator of the evolutionary dynamics of asexual populations with explicit genomic sequences. SodaPop is written in an object-oriented programming (OOP) framework, where the effects of population structure and the biophysical effects of mutations can be explored simultaneously. The input is the spectrum of mutational effects on biochemical or biophysical properties, which can be derived from protein engineering methods (Jia et al., 2015; Kumar et al., 2006; Laimer et al., 2015; Yin et al., 2007) or from exhaustive mutagenesis experiments such as deep mutational scanning (DMS) (Bloom, 2014; Firnberg et al., 2014; Fowler and Fields, 2014). SodaPop also allows flexibility in defining fitness functions from biochemical/biophysical models that describe the evolution of proteins. Finally, the OOP framework facilitates the integration of new attributes and parameters into the model, as well as the customization of input and output data.
To the best of our knowledge, SodaPop is the first publicly available tool that explicitly combines the role of protein biophysics and population dynamics. The main program is implemented in C++ with a command line interface. We also provide scripts for analysis and visualization of the simulation results. Source code, binaries and documentation can be downloaded freely from https://github.com/louisgt/SodaPop under the GNU GPLv3 license. This software is portable on any POSIX-compliant operating system, including Linux and Mac OS/X, or on Windows using the Cygwin environment.
2 Materials and methods
2.1 Hierarchical and object-oriented description of an asexual population
We designed a hierarchical and object-oriented representation of an asexual population (Fig. 1A). A population of N cells is contained in a dynamic array. Cell is a data structure defined by the Cell.h class whose attributes include its fitness (the reproductive capacity and defined in greater detail below), its mutation rate, and the array of genes in the genome (Supplementary Fig. S1). Each cell is also assigned an ID (15-nt barcode), which can be used to track different lineages in the population. The Cell data structure can be extended to include other attributes and functions. Additionally, each Cell contains a vector of Genes, a data structure defined by the Gene.h class (Supplementary Fig. S2). Each Gene is characterized by its nucleotide sequence, thermodynamic folding stability, enzymatic efficiency, intracellular abundance and essentiality. Each instance of a gene also contains the variables Na and Ns, which track the non-synonymous and synonymous substitutions, respectively. The Gene.h class is also extendable to include other features such as protein-protein interaction and information on binding partners (see Supplementary Manual, Section 8.2).
2.2 Evolutionary algorithm
SodaPop uses an adapted Wright-Fisher model with selection (Fisher, 1922; Wright, 1931) to evolve the population in discrete and non-overlapping generations (Fig. 1C). To propagate to the next generation, each cell (the parent) in the population is replicated to k progenies. The number k is drawn from a binomial distribution with N trials and a probability of success equal to the relative fitness of the parent wc. This value is defined as the fitness of the parent cell normalized by the sum of the fitness of all cells in the population (Fig. 1C). The total number of progenies N could differ from the target population size Ne, thus to maintain a constant population size, it is scaled appropriately between generations. If N < Ne, randomly sampled cells from the offsprings are replicated. If N > Ne, ΔN randomly chosen progenies are removed. Population size can also be adjusted between generations to mimic population bottleneck experiments or expansion into an ecological niche (Barrick and Lenski, 2013; Ebert, 1998; Gullberg et al., 2011).
During replication, the genome of each new progeny can acquire mutations at a rate Lμ, where μ is the mutation rate per base pair per generation and L is the length of the genome. Specifically, for each replicated cell, the number of mutations m is drawn from a binomial distribution with L trials and a probability of success equal to μ. All mutations have equal likelihood of occurring anywhere in the genome. Each mutation is thus randomly mapped to a specific site in a particular gene. Depending on the fitness landscape chosen or defined by the user (Fig. 1B), the effect of a mutation on the biophysical properties of the gene product and fitness can either come from experiments such as deep mutational scan, from physical force field calculations, or from a predefined distribution (Supplementary Manual, Chapter 5). The mutation model is also reversible, so that all effects on protein properties contributed by the previous allele are removed, to be replaced by those of the incoming substitution. The fitness landscape models are described in Section 2.3.
In the course of the simulation, SodaPop saves a snapshot of the population every T generations, as defined by the user. The snapshot of all the genomic sequences of cells in the population allows for the reconstruction of evolutionary trajectories and the calculation of evolutionary rates. Additionally, this snapshot can be used as input for subsequent simulations. Users can also choose from multiple output formats according to their specific needs. For instance, a shorter output format only includes information at the level of cells, namely, barcodes and fitness values. A more complete format includes gene information for every cell, including DNA and amino acid sequences. This functionality allows users to tune both the level of detail required and the speed and memory usage of the program. For complete documentation on output formats, refer to the Supplementary Manual Section 3.4.
SodaPop is built upon memory-efficient data structures and a fast algorithm to achieve high computational performance and to minimize the general trade-off between flexibility and runtime (Carvajal-Rodriguez, 2008). The program can readily execute simulations in the order of 106 individually defined cells with runtimes clocking under a few hours on an ordinary four-core desktop computer with 8 or 16GB RAM.
2.3 Fitness landscapes
SodaPop allows users to choose from several fitness landscapes based on protein biochemistry. These fitness landscapes have been used to model and explain the rates of protein evolution (Bloom et al., 2007; Serohijos et al., 2012; Taverna and Goldstein, 2002), polymorphisms in protein coding regions (Serohijos and Shakhnovich, 2014a, b), epistasis (Bershtein et al., 2006; Bloom et al., 2007) and the log-normal distribution of protein evolutionary rates in genomes (Wolf et al., 2009).
2.4 Fitness effects of mutations
When a random mutation occurs, it may change the biochemical and biophysical properties of a protein and, in turn, the fitness of the cell. SodaPop has three approaches to model the fitness effects of mutations (Fig. 1b).
2.4.1 Mutational effects derived from a distribution
The effects of random mutations on the folding stability of globular proteins can be characterized as Gaussian distribution with mean μ = 0.6 kcal/mol and standard deviation σ = 0.9 kcal/mol (Tokuriki et al., 2007). These results arose from both comprehensive mutagenesis of several proteins and then estimating the effects of mutations on stability using physical force fields (Tokuriki et al., 2007). These distributions are also in agreement with >5000 stability measurements of purified proteins (Kumar et al., 2006). Thus, the user can define the DFE as a two-parameter distribution of the form N(μ, σ) or in the command-line (Supplementary Manual, Section 5.1). Instead of drawing ΔΔG values from a Gaussian, a Gamma distribution can be used to draw selection coefficients (Eyre-Walker et al., 2006; Nielsen and Yang, 2003; Tamuri et al., 2012) and calculate the corresponding fitness.
2.4.2 Mutational effects derived from physical force fields
Users can also provide as input a matrix that describes the change in folding stability () or activity () for all possible one-away substitutions at all residues. These quantities can be provided by the user (using parameter –i) as look-up tables that are accessed at runtime. The entries in these tables are derived from computational protein engineering tools such as Rosetta (Das and Baker, 2008), Eris (Yin et al., 2007) or FoldX (Guerois et al., 2002) prior to the evolutionary simulation. Updating the protein folding stability or activity also updates the fitness using either Equations (2)–(4).
2.4.3 Mutational effects derived from experiment
The user can also provide as input fitness effects from experimental DMS approaches (Araya and Fowler, 2011). DMS combines mutational library generation, selection, and high-throughput sequencing to assay the fitness of up to 95% of possible one-away mutations to a protein. To date, over a dozen systematic and exhaustive DMS assays have been conducted in various proteins to determine their local fitness landscape (Fowler and Fields, 2014; Wrenbeck et al., 2017a). New computational tools can leverage the information in these experimental datasets to predict mutational effects for full proteomes (Gray et al., 2018), allowing the creation of comprehensive substitution matrices for thousands of proteins. Users can input these matrices using the parameter −i (Supplementary Manual, Section 3.3).
Indels are another major source of innovation in protein evolution. Although the current version of SodaPop does not explicitly handle them, some types of indels can already be modeled (see Section 8.3 of the Supplementary Manual).
2.5 Analysis of simulation data
The SodaPop package includes shell and R scripts for visualizing the population dynamics and the time-series of the fitness and biophysical properties of proteins. These scripts also allow for the tracking of lineages and clonal structure and the calculation of protein evolutionary rates. A detailed list and description of outputs and analyses can be found in Chapter 6 of the User Manual. Since the program registers explicit DNA and protein sequences for all the cells in the population, a single simulation can generate thousands to millions of sequences per saved generation. These sequences can be used as starting points for high-resolution methods in molecular evolution and phylogenetics. For example, the sequences can be used to create multiple sequence alignments and apply standard molecular evolution analyses, such as the McDonald-Kreitman test or Tajima’s D.
3 Results
We describe four applications of SodaPop to demonstrate its broad scope and flexibility.
3.1 Application I: population dynamics on fitness landscapes based on protein biochemistry
To show how SodaPop can be used to explore population dynamics and protein properties on biochemical fitness landscapes, we simulated a population of Ne = 104 cells on the fitness landscape defined by Equation (3). Each cell contains ten genes in the folate biosynthesis pathway of Escherichia coli. Shown in Figure 2A and B are the mean population fitness and the average stability of the ten genes over the course of the simulation. At a finer temporal resolution, we can trace the segregation and eventual fixation of arising mutations and the effects of clonal interference (Fig. 2C;Supplementary Fig. S4).
3.2 Application II: simulating population dynamics with barcoded cells for lineage-tracking and competition assays
The ability to track the lineages of cells and frequency of clones in an evolving population is crucial because it enables estimating the selective advantage of mutations, the extent of clonal interference and the establishment time of adaptive mutations. Tracking the population dynamics typically makes use of neutral genetic or fluorescent markers (Hegreness et al., 2006; Illingworth and Mustonen, 2012; Moura de Sousa et al., 2013; Pinkel, 2007; Zhang et al., 2012). More recently, through the introduction of randomized and unique barcodes into the chromosomes, it is now possible to track the population dynamics and lineages at the resolution of single cells [Blundell and Levy (2014); Gerrits et al. (2010); Levy et al. (2015); Venkataram et al. (2016)]. Nonetheless, because the resulting dynamics can be complex, these experimental results are complemented by simulations that could provide null expectations. SodaPop allows for this functionality. By assigning a unique barcode identifier to each cell, one can trace the lineage and history of cells in the simulation. As an example, we performed a simulation of laboratory evolution of Ne = 104 cells, each assigned a 15-nt barcode (Fig. 3). To mimic standing genetic variation, different fitness values can be assigned to the barcodes (Supplementary Manual Section 4.2). Figure 3A is a lineage density plot showing the relative share of each barcode in the population and Figure 3B shows the frequency of each barcode.
3.3 Application III: simulating coding sequence evolution under the constraints of both population dynamics and selection for folding stability
Next, we demonstrate how SodaPop can be used to perform simulations of protein sequence evolution that account for both effects of selection for biophysical properties, such as folding stability, and population dynamics. Specifically, we show the performance of SodaPop in recapitulating the amino acid conservation among orthologs. As a model system, we chose aminodeoxychorismate synthase (pabB gene in Escherichia coli), an enzyme in the folate biosynthesis pathway.
To generate simulated orthologs, we performed evolutionary simulations of 105 cells with E. coli pabB for 106 generations, under selection for folding stability (Eq. 3). The effects of random mutations on fitness are derived from a 452 residue × 20 amino acid matrix that contains values for changes in folding stability, . The values are estimated using a physical force field (Yin et al., 2007) and using the pabB 3D structure from Protein Data Bank (PDB ‘1K0E’). We first performed an equilibration simulation where an initially monoclonal population was evolved for 105 generations to reach mutation-selection balance. Then, to mimic divergence from a common ancestor, we used the endpoint of the equilibration simulation as the starting population for 250 independent evolutionary simulations. These divergent simulations are performed for 106 generations, thereby ensuring that the distribution of pairwise sequence identities for simulated sequences matches that of the extant orthologs of pabB.
To compare the simulated sequences with extant orthologs of pabB, we retrieved the top hits of a protein sequence search in OrthoDB for bacteria (Waterhouse et al., 2013). We excluded sequences with an absolute length difference of more than 20 bp. We then aligned the 657 remaining sequences with hmmalign (http://hmmer.org/) to the corresponding HMM of Pfam domains in pabB (PF04715 and PF00425) and trimmed the flanking gaps from the alignment. The amino acid conservation of simulated and orthologous pabB sequences from multiple sequence alignment is shown in Supplementary Figures S5 and S6, respectively.
Next, we explored if simulations with SodaPop could also capture the co-evolution between sites and whether this result is dependent on Ne. Here, co-evolution is calculated by mutual information (MI) using MISTIC (Simonetti et al., 2013). Indeed, increasing the stringency of selection due to larger Ne increases the number of co-evolving sites and the strength co-evolution between those sites (Fig. 5A). Interestingly, population dynamics seems to influence the rank-ordered distribution of cumulative site co-evolution (CMI) with increasing Ne (Fig. 5B). In protein engineering, site co-evolution is now being used to predict de novo 3D structure of proteins and protein complexes (Marks et al., 2012). Being able to explore how this property is influenced by population dynamics will be of practical importance to the community.
3.4 Application IV: simulating sequence evolution with selection for both folding stability and catalytic activity
Mutations affect not only protein stability, but also enzyme activity. For some residues and structural motifs, the trade-off between activity and stability can constrain evolutionary paths (Meiering et al., 1992). Thus, we compared two models of enzyme evolution with SodaPop: (1) a model with effects on protein stability alone and (2) a model with effects on enzyme activity in addition to stability. For the first model, we used FoldX’s PositionScan command (Guerois et al., 2002) to get the landscape of effects on folding stability of a protein. For the second model, we took advantage of a recent deep-mutational scanning assay comprising >96.3% of all possible one-away non-synonymous substitutions for hydrolysis activity of the amiE amidase (Wrenbeck et al., 2017b). We regressed the previously calculated effects on folding from the DMS data to extract the contribution of each mutation to enzyme activity (Supplementary Methods). We find that supplementing the effects on protein stability with effects on the catalytic activity increases the recapitulation of sequence conservation in natural orthologs in comparison with sequences from simulation based on folding stability alone (Table 1, Supplementary Fig. S7).
. | Pearson correlation . | P-value . |
---|---|---|
Selection for folding | 0.33 | 3.85 × 10-10 |
Selection for folding and activity | 0.43 | 2.2 × 10-16 |
. | Pearson correlation . | P-value . |
---|---|---|
Selection for folding | 0.33 | 3.85 × 10-10 |
Selection for folding and activity | 0.43 | 2.2 × 10-16 |
. | Pearson correlation . | P-value . |
---|---|---|
Selection for folding | 0.33 | 3.85 × 10-10 |
Selection for folding and activity | 0.43 | 2.2 × 10-16 |
. | Pearson correlation . | P-value . |
---|---|---|
Selection for folding | 0.33 | 3.85 × 10-10 |
Selection for folding and activity | 0.43 | 2.2 × 10-16 |
3.5 Performance and runtime
We benchmarked SodaPop for six population sizes spanning five orders of magnitude with cells containing ten genes (total genome size 6.5 kbp). All simulations were run on a standard iMac desktop with a 3.2GHz Intel Core i5 processor and 16GB RAM. Figure 6 shows that runtime is quasi-monomial with respect to population size, with an exponent of 1.1 (Supplementary Fig. S8). Simulating up to a million cells for long time periods is entirely tractable using standard desktop computers. We limited our desktop benchmarking to N = 106 cells, as higher orders of magnitude introduce a shift in performance due to a RAM bottleneck. The simulation of populations with higher orders of magnitude requires a larger amount of memory than the current standard in commercial desktop computers. However, these larger population sizes can be simulated on high-performance computing clusters where memory allocation is not limiting.
4 Discussion and conclusion
We have created SodaPop, a fast and scalable tool for evolutionary simulations based on biochemical and biophysical fitness landscapes. Considering the need to address questions at the interface of molecular evolution and population genetics, and with most of the current computational methods unable to account for explicit clonal dynamics, we believe SodaPop provides a comprehensive and extensible framework that can encompass a wide array of evolutionary scenarios. Briefly, we showed that our program can be used, among other things, to simulate evolution under selection for specific biophysical properties, to provide null expectations for lineage-tracking experiments of laboratory evolution, and to model coding sequence evolution and co-evolution. To broaden these possibilities and achieve greater predictive power, we intend to implement future models of fitness landscapes, such as a systems-level description of the cell that could integrate omics data—transcription level, gene regulation networks and protein-protein interactions. We also aim to expand the scope of the program to other types of genetic changes by integrating features of recombination and indels to SodaPop. In the case of recombination, the first and simplest approach is to assume an ad hoc distribution of effects, as done in other forward simulation algorithms. This approach is appropriate for modeling whole genome recombination events (inter-genic). However, this does not account for the consequence of recombination in specific genes (intra-genic) and on their biophysical properties, which is the overaching motivation for SodaPop. The second approach is to phenomenologically model the fitness effect of recombination and its dependence on recombination sites within proteins. Random recombination of homologous sequences found that crossover sites in the middle of the protein are more deleterious because it likely impacts the packing of the core (Romero and Arnold, 2009). The third approach is to have a sequence-based model of recombination and to calculate the biophysical properties of the recombinant sequences. Considering that for a given evolutionary run, at least 106 sequences need to be folded (but few of which will be selected), on-the-fly calculation of folding stability for a sequence is possible only in lattice models of protein folding (Shakhnovich and Gutin, 1993; Voigt et al., 2002), barring the associated computational cost.
In the case of indels, their effects in protein evolution are more nuanced and depend on structural and biophysical considerations. For example, random indels occurring in reverse turns, loops, or surfaces are considered to be less deleterious than those in beta sheets, helices, or protein cores (Benner et al., 1993; Hsing and Cherkasov, 2008; Pascarella and Argos, 1992). A detailed model of indels will require explicit folding simulation of proteins for each mutation, which can also be performed using lattice models of proteins.
The object-oriented design of SodaPop will facilitate these future developments. In a technical perspective, there are several features that could improve the performance and practical use of SodaPop. First, modeling mutation history for each cell using a linked list of mutations instead of explicit sequence change could reduce the memory required to store cell objects by a significant factor without incurring any information loss. Second, circumventing the command-line application with a graphical user interface (GUI) wrapper will facilitate user interaction in creating input files, setting up simulations and choosing the appropriate parameters. Third, implementing multi-threading options in the main program loop and in the analysis pipeline will allow the program to run on multiple processors in parallel, which can significantly improve runtime on high-performance computing clusters. These features are currently under development for future versions.
Acknowledgements
We would like to thank members of the Serohijos Lab for testing the program and contributing valuable questions and ideas.
Funding
A.W.R.S. acknowledges a grant from (Natural Sciences and Engineering Research Council) and start-up funds from Université de Montréal. L.G. is a recipient of a doctoral fellowship from Université de Montréal's Faculté des études supérieures et postdoctorales. R.D.F. was supported by ENSEIRB-MATMECA Bordeaux.
Conflict of Interest: none declared.
References