Amir Yassin
email: amir.yassin@universite-paris-saclay.fr

Nelly Gidaszewski

Vincent Debat

Jean R David

Long-term evolution of quantitative traits in the Drosophila melanogaster species subgroup

Keywords: phylogenetic constraints, modularity, pigmentation, body size, bristle number, ovariole number

Quantitative genetics aims at untangling the genetic and environmental effects on phenotypic variation. Trait heritability, which summarizes the relative importance of genetic effects, is estimated at the intraspecific level, but theory predicts that heritability could influence long-term evolution of quantitative traits. The phylogenetic signal concept bears resemblance to heritability and it has often been called species-level heritability. Under certain conditions, such as trait neutrality or contribution to phylogenesis, within-species heritability and between-species phylogenetic signal should be correlated. Here, we investigate the potential relationship between these two concepts by examining the evolution of multiple morphological traits for which heritability has been estimated in Drosophila melanogaster. Specifically, we analysed 42 morphological traits in both sexes on a phylogeny inferred from 22 nuclear genes for nine species of the melanogaster subgroup. We used Pagel's λ as a measurement of phylogenetic signal because it is the least influenced by the number of analysed taxa. Pigmentation traits showed the strongest concordance with the phylogeny, but no correlation was found between phylogenetic signal and heritability estimates mined from the literature. We obtained data for multiple climatic variables inferred from the geographical distribution of each species. Phylogenetic regression of quantitative traits on climatic variables showed a significantly positive correlation with heritability. Convergent selection, the response to which depends on the trait heritability, may have led to the null association between phylogenetic signal and heritability for morphological traits in Drosophila. We discuss the possible causes of discrepancy between both statistics and caution against their confusion in evolutionary biology.

Introduction

Heritability is a central concept in genetics with crucial applications in medicine and agriculture [START_REF] Visscher | Heritability in the genomics era -concepts and misconceptions[END_REF]. It refers to the tendency of genetically related individuals (e.g., parent-offsprings, sibs, etc.) to resemble each other for a particular phenotype more than they do with other individuals randomly drawn from the population, regardless of environmental conditions.

The heritability estimate for a trait is a strong indicator for the response of that trait to selection [START_REF] Falconer | Introduction to Quantitative Genetics[END_REF][START_REF] Lynch | Mutations and quantitative genetic variation: lessons from Drosophila[END_REF]. Consequently, long-term parallel selective regimes in different populations could lead to similarity in phenotypic values for a heritable trait among genetically-distant individuals, i.e. convergent evolution, due to selection on either independent or shared alleles (e.g., [START_REF] Gross | A Novel Role for Mc1r in the Parallel Evolution of Depigmentation in Independent Populations of the Cavefish Astyanax mexicanus[END_REF][START_REF] Jones | The genomic basis of adaptive evolution in threespine sticklebacks[END_REF][START_REF] Ralph | Convergent Evolution During Local Adaptation to Patchy Landscapes[END_REF]. Therefore, heritability is usually measured within populations using a limited number of generations, where relatedness could be traced and experimental conditions could be controlled.

In spite of such experimental limitations, the concept of heritability has been extended to comparative approaches studying trait evolution across species, where it has mostly been measured in term of phylogenetic signal and was called species-level or phylogenetic heritability (e.g., [START_REF] Jablonski | Heritability at the Species Level: Analysis of Geographic Ranges of Cretaceous Mollusks[END_REF][START_REF] Lynch | Methods for the Analysis of Comparative Data in Evolutionary Biology[END_REF][START_REF] Blomberg | Testing for Phylogenetic Signal in Comparative Data: Behavioral Traits Are More Labile[END_REF][START_REF] Queiroz A De | The Phylogeny of a Species-Level Tendency: Species Heritability and Possible Deep Origins of Bergmann's Rule in Tetrapods[END_REF][START_REF] Waldron | Null Models of Geographic Range Size Evolution Reaffirm Its Heritability[END_REF][START_REF] Cardini | Does the skull carry a phylogenetic signal? Evolution and modularity in the guenons[END_REF][START_REF] Shirreff | How effectively can HIV phylogenies be used to measure heritability?[END_REF][START_REF] De Villemereuil | General Quantitative Genetic Methods for Comparative Biology[END_REF][START_REF] Vrancken | Simultaneously estimating evolutionary history and repeated traits phylogenetic signal: applications to viral and host phenotypic evolution[END_REF][START_REF] Houle | Mutation predicts 40 million years of fly wing evolution[END_REF][START_REF] Oliver | Phylogenetic signal in tooth wear dietary niche proxies[END_REF][START_REF] Scher | Application of remote sensing technology to estimate productivity and assess phylogenetic heritability[END_REF]. In phylogenetics, a trait with a phylogenetic signal is a trait showing a tendency of closely-related species to resemble each other more than to other species randomly drawn from the phylogeny [START_REF] Pagel | Detecting Correlated Evolution on Phylogenies: A General Method for the Comparative Analysis of Discrete Characters[END_REF][START_REF] Revell | Phylogenetic Signal, Evolutionary Process, and Rate[END_REF][START_REF] Münkemüller | How to measure and test phylogenetic signal[END_REF], a definition that closely resembles that of heritability given above, species playing the role of individuals. However, unlike heritability, phylogenetic signal is by definition measured in situations where convergent selection or genetic introgression could not be controlled, and its biological interpretation has long been problematic [START_REF] Losos | Seeing the Forest for the Trees: The Limitations of Phylogenies in Comparative Biology[END_REF]. Consequently, the relationship between short-term quantitative genetics estimates of heritability and long-term phylogenetic signals remains unclear. Indeed, a positive correlation is expected for neutrally evolving traits, but the relationship is more complex for traits related to fitness.

For nearly a century, Drosophila melanogaster has been a quantitative genetics paradigm [START_REF] Roff | Quantitative genetics and fitness: lessons from Drosophila[END_REF]Mackay 2010). At the morphological level, four trait categories, namely size, shape, setation (bristle number) and pigmentation have been thoroughly studied on several anatomical structures (e.g., wing, thorax, different abdominal segments, etc.). This has resulted into a wealth of data on heritability, modularity, sexual dimorphism, phenotypic plasticity, response to selection and geographical variation of those traits (reviewed in [START_REF] Gibert | Comparative analysis of morphological traits among Drosophila melanogaster and D. simulans: genetic variability, clines and phenotypic plasticity. In: Drosophila melanogaster, Drosophila simulans: So Similar, So Different[END_REF]. With the advent of genome-wide association techniques, the genomic basis of many of traits is now being unravelled [START_REF] Mackay | Drosophila bristles and the nature of quantitative genetic variation[END_REF][START_REF] Turner | Population-Based Resequencing of Experimentally Evolved Populations Reveals the Genetic Basis of Body Size Variation in Drosophila melanogaster[END_REF][START_REF] Takahashi | Deficiency Screening for Genomic Regions with Effects on Environmental Sensitivity of the Sensory Bristles of Drosophila Melanogaster[END_REF][START_REF] Bastide | A genome-wide, fine-scale map of natural pigmentation variation in Drosophila melanogaster[END_REF][START_REF] Bastide | A Variable Genetic Architecture of Melanic Evolution in Drosophila melanogaster[END_REF][START_REF] Dembeck | Genetic Architecture of Abdominal Pigmentation in Drosophila melanogaster[END_REF][START_REF] Vonesch | Genome-Wide Analysis Reveals Novel Regulators of Growth in Drosophila melanogaster[END_REF][START_REF] Lobell | The Genetic Architecture of Ovariole Number in Drosophila melanogaster: Genes with Major, Quantitative, and Pleiotropic Effects[END_REF][START_REF] Lafuente | Genetic basis of thermal plasticity variation in Drosophila melanogaster body size[END_REF][START_REF] Pitchers | A Multivariate Genome-Wide Association Study of Wing Shape in Drosophila melanogaster[END_REF]. Of these traits, only wing shape has been analysed using geometrical multivariate approaches within D. melanogaster (e.g., [START_REF] Debat | Allometric and Nonallometric Components of Drosophila Wing Shape Respond Differently to Developmental Temperature[END_REF][START_REF] Pitchers | Altitudinal Clinal Variation in Wing Size and Shape in African Drosophila Melanogaster: One Cline or Many?[END_REF] as well as among the nine species of the melanogaster subgroup [START_REF] Gidaszewski | Evolution of sexual dimorphism of wing shape in the Drosophila melanogaster subgroup[END_REF][START_REF] Klingenberg | Testing and quantifying phylogenetic signals and homoplasy in morphometric data[END_REF][START_REF] Houle | Mutation predicts 40 million years of fly wing evolution[END_REF]. However, the relationship between short-term heritability and phylogenetic signals in wing shape traits is still unknown, perhaps due to the nature of this character which could convolute homologization of multivariate statistical summaries (e.g., principal components) within and between species.

In the present study, we investigate the evolution of 42 morphological traits related to size, shape, bristle number, pigmentation and ovariole number, in the nine species of the Drosophila melanogaster subgroup. Heritability of most of those traits has been estimated in D. melanogaster in terms of intraclass correlation (t) using the same isofemale line approach [START_REF] Delpuech | Phenotypic Plasticity and Reaction Norms in Temperate and Tropical Populations of Drosophila Melanogaster: Ovarian Size and Developmental Temperature[END_REF][START_REF] Gibert | Describing the Evolution of Reaction Norm Shape: Body Pigmentation in Drosophila[END_REF][START_REF] Gibert | Phenotypic plasticity of abdomen pigmentation in two geographic populations of Drosophila melanogaster: male-female comparison and sexual dimorphism[END_REF][START_REF] Karan | Body size reaction norms in Drosophila melanogaster: temporal stability and genetic architecture in a natural population[END_REF][START_REF] Araripe | Divergent abdominal bristle patterns in two distantly related drosophilids: antero-posterior variations and sexual dimorphism in a modular trait[END_REF]). Although t involves in theory nonadditive components, the ratio between t and independent h 2 estimates of heritability was found in practice to be 0.95, indicating that t is strongly related to h 2 [START_REF] David | Isofemale lines in Drosophila: an empirical approach to quantitative trait analysis in natural populations[END_REF]. Moreover, t estimates for body size and bristle number traits in the distantly-related drosophilid genus Zaprionus approaches those inferred for D. melanogaster [START_REF] David | Quantitative trait analysis and geographic variability of natural populations of Zaprionus indianus, a recent invader in Brazil[END_REF][START_REF] Araripe | Divergent abdominal bristle patterns in two distantly related drosophilids: antero-posterior variations and sexual dimorphism in a modular trait[END_REF], suggesting that heritability estimates in D. melanogaster are likely a good indicator for the trait heritability among its closely-related species. We inferred a molecular phylogeny for the nine species using 22 nuclear genes, and estimated the phylogenetic signal of each trait using Pagel's λ [START_REF] Pagel | Detecting Correlated Evolution on Phylogenies: A General Method for the Comparative Analysis of Discrete Characters[END_REF], because this statistic is the least influenced by the number of analysed species [START_REF] Münkemüller | How to measure and test phylogenetic signal[END_REF]). We also obtained data for 14 climatic variables from the geographical distribution of each species. We found that t did not correlate with λ, but it correlated with the phylogenetic regression on the climatic variables. Because heritability determines response to selection, our results indicate that convergent evolution has likely shaped long-term evolution of classical quantitative traits in Drosophila leading to the dissociation between their within-species heritability estimates and their between-species phylogenetic signal.

Materials and Methods

Drosophila strains

With two exceptions, indicated below, we always used mass cultures founded by a set of at least 30 females or isofemale lines. For the two cosmopolitan species, D. melanogaster and D. Groups of 10-20 adult flies were isolated on usual cornmeal sugar food and kept there for 3-4 days. They were then introduced into rearing vials containing a high nutrient, killed yeast medium for producing the experimental flies at a temperature of 21°C. Each group of adults was transferred to a fresh vial for at least 3 successive days. With this procedure, population density in a vial was less than 100. After emergence, young adults were transferred to fresh vials and kept at least for 3 days before being measured. Measurements were made with a binocular microscope equipped with an ocular micrometer, on anesthetized flies (n = 10 per sex per strain). For anesthesia we used triethylamine (also called flynap), which has the advantage of producing a long lasting anesthesia on immobile flies.

Traits measured

A total of 20 different traits were measured in females and 15 in males. They correspond either to linear or area measures (metric traits) or to numbers (meristic traits). Also 4 morphological indices or ratios, were considered in both sexes. These traits were measured either on the thorax or on the abdomen.

Body size traits. Three linear measurements were taken, and the dimensions transformed into 100x mm: wing length (WL), from the thoracic articulation to the tip of the wing; thorax length (TL) from the neck to the tip of the scutellum, on a left side lateral view; and thorax width (TW), as the distance between the two posterior large sternopleural bristles, from a ventral view.

Body shape indices.

From the above data we calculated, on each fly, three ratios: WL/TL which is inversely proportional to the wing loading [START_REF] Pétavy | Growth temperature and phenotypic plasticity in two Drosophila sibling species: probable adaptive changes in flight capacities[END_REF]; TL/TW which is the thorax elongation index; and WL/TW which also provides information on wing elongation.

Bristles number. Bristles were counted on the thorax and on the ventral side of the abdomen.

Thoracic, sternopleural bristles were counted on each side and the sum (STPn) was used in calculations. We also measured the length of the two major bristles, anterior (STPa) and posterior (STPp) sternopleurals, and also calculated their ratio (STPp/a). The inverse of this ratio, called spernopleural index [START_REF] Sturtevant | The classification of the genus Drosophila, with descriptions of nine new species[END_REF], is often used in systematics. Abdominal bristles were counted on the successive sternites, from segment 2 up to segment 5 in males (A2-A5) and segment 7 in females (A2-A7). The bristles were also counted on each side of tergite 7 in females, and the sum calculated (T7).

Abdominal pigmentation. This was measured on successive tergites (P2 to P7 in females and P2 to P6 in males). The surface of black area at the posterior part of each tergite was visually scored by establishing 11 phenotypic classes from 0 (no black pigment) up to 10 (tergite completely black) [START_REF] David | Abdominal pigmentation and growth temperature in Drosophila melanogaster: Similarities and differences in the norms of reaction of successive segments[END_REF]).

Ovariole number. Mature females were dissected and the ovariole number counted on each ovary and then the sum calculated (OV).

Phylogeny inference

We chose the 22 nuclear genes that were analysed by [START_REF] Akashi | Molecular Evolution in the Drosophila melanogaster Species Subgroup: Frequent Parameter Fluctuations on the Timescale of Molecular Divergence[END_REF] to investigate codon usage evolution in five species of the melanogaster subgroup: D. melanogaster, D. simulans, D. teissieri, D. erecta and D. orena. We verified sequences of these genes in the five species of the subgroup with published, annotated complete genome sequence: D. melanogaster, D. simulans, D. sechellia, D. yakuba and D. erecta [START_REF] Clark | Evolution of genes and genomes on the Drosophila phylogeny[END_REF]). We also obtained their DNA sequences from the published genomes of D. mauritiana [START_REF] Nolte | Genome-wide patterns of natural variation reveal strong selective sweeps and ongoing genomic conflict in Drosophila mauritiana[END_REF]) and D. santomea [START_REF] Turissini | The evolution of reproductive isolation in the Drosophila yakuba complex of species[END_REF] using BLAST [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF]. Sequences from each gene were separately aligned with Muscle [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF]) as implemented in the MEGA 7.0 software package [START_REF] Kumar | MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets[END_REF]. We also used MEGA to estimate the best substitution model for each gene, which were HKY + Γ for Adh, Adhr, Amyrel, boss,dpp,Fur2,g,Gld,Gpdh,Mcm3,pav,per,RpII215,run,twi and Zw, HKY + I for Amy-p and Cyp28c1, GTR + Γ for LanB2 and ry, and HKY + Γ + I for Osbp. We concatenated the genes (41,382 nucleotides) and inferred a Bayesian phylogeny using MrBayes v3.2.4 [START_REF] Ronquist | MrBayes 3.2: Efficient Bayesian Phylogenetic Inference and Model Choice Across a Large Model Space[END_REF], applying the estimated substitution model to each gene partition. Two runs of 500,000 generations were conducted and sampled every 1,000 generations under a strict clock model. We assessed convergence -the average standard deviation of split frequency < 0.01 and the potential scale reduction factor ~1.00 -using MrBayes. A burn-in period of 25% of samples was used. Alignment and MrBayes commands are given in Supplementary Text 1.

Phylogenetic signal

Morphometric data were analysed as univariate variables to infer the mean and standard error per species in R (www.r-project.org). R was also used to conduct Principal Component Analysis (PCA) on the complete (multivariate) data set after excluding ratio traits. Mean values were mapped on the species phylogeny (excluding data of temperate strains of D. melanogaster and D. simulans). We inferred phylogenetic signal (species-level heritability) in terms of [START_REF] Pagel | Detecting Correlated Evolution on Phylogenies: A General Method for the Comparative Analysis of Discrete Characters[END_REF] lambda (λ). We used the Generalized Least Square (GLS) approach to infer λ as implemented in the Continuous Model A method assuming random walk in the BayesTraits v2.0 software package [START_REF] Pagel | BayesTraits V2.0. Computer program and documentation[END_REF]. All analyses were conducted using MCMC chains of 1,010,000 iterations and 1,000 post-burnin trees on trees generated by MrBayes to account for topological uncertainties.

The harmonic mean of the likelihoods, i.e., marginal likelihoods of the MCMC runs, was checked for stability. We estimated λ significance in terms of Bayes factors (BFs) by comparisons with runs where λ was restricted to 0 (i.e. the trait evolution is independent from the phylogeny) or 1 (i.e. the phylogeny correctly predicts the trait evolution). BF was estimated as 2 * (marginal likelihood of the best-fit model -marginal likelihood of the worse-fit model). We considered BF > 2, 5 or 10 as supportive, strong, or very strong support of the best-fit model, respectively, following the general convention [START_REF] Pagel | BayesTraits V2.0. Computer program and documentation[END_REF].

In order to examine the relationship between phylogenetic signal (λ) and short-term heritability, we mined the quantitative genetics literature in D. melanogaster to obtain heritability estimates for the morphometrical traits. Because heritability could be estimated using various approaches, we opted for the sake of comparison for a single methodology, i.e. intra-class correlation (t), that has been used in most of the traits analysed here [START_REF] Delpuech | Phenotypic Plasticity and Reaction Norms in Temperate and Tropical Populations of Drosophila Melanogaster: Ovarian Size and Developmental Temperature[END_REF][START_REF] Gibert | Describing the Evolution of Reaction Norm Shape: Body Pigmentation in Drosophila[END_REF][START_REF] Gibert | Phenotypic plasticity of abdomen pigmentation in two geographic populations of Drosophila melanogaster: male-female comparison and sexual dimorphism[END_REF][START_REF] Karan | Body size reaction norms in Drosophila melanogaster: temporal stability and genetic architecture in a natural population[END_REF][START_REF] Araripe | Divergent abdominal bristle patterns in two distantly related drosophilids: antero-posterior variations and sexual dimorphism in a modular trait[END_REF]. To the best of our knowledge, no heritability estimates using t was made for any of the analysed traits in other species of the melanogaster subgroup. However, t estimates for body size and bristle number traits in the distantly-related drosophilid genus Zaprionus approaches those inferred for D. melanogaster [START_REF] David | Quantitative trait analysis and geographic variability of natural populations of Zaprionus indianus, a recent invader in Brazil[END_REF][START_REF] Araripe | Divergent abdominal bristle patterns in two distantly related drosophilids: antero-posterior variations and sexual dimorphism in a modular trait[END_REF], suggesting that heritability estimates in D. melanogaster are likely a good indicator for the trait heritability among its closely-related species. Therefore, we inferred the coefficient of correlation between λ and t, as an indication for the relationship between short-term heritability and long-term phylogenetic signal.

Phylogenetic regression on climatic variables

In order to determine which environmental variables might also explain the evolution of each trait, we conducted a regression analysis for each trait (except for shape ratios) on species means of 14 climatic variables. These included seven climatic sets (namely, average (Tavg), minimum (Tmin) and maximum (Tmax) air temperatures and earth temperature (TEarth) measured in °C, relative humidity (RH) in %, atmospheric pressure (AP) in kPa, and wind speed (WS) at 50 m above ground in m/s) and seven radiation sets (namely, insolation incident (II), diffuse radiation (DR), direct normal radiation (DNR) and latitude tilt radiation (LTR), each measured in kWh/m 2 /day, clear sky (CS) in number of days with no clouds, no sun (NS) in number of days with clouds, and ultra-violet index (UVI) in units of 25 mW/m 2 ). We first obtained geographical coordinates (i.e. latitudes and longitudes) from TaxoDros v1.04 (Bächli 2016; last accessed December 2020) for each locality from which each of the nine species was collected in tropical Africa (i.e. south to latitude 17°N). We then extracted for each locality climatic data averaged over 22 years (from 1983 to 2005) as monthly averages from the NASA Surface meteorology and Solar Energy: Global Data Sets website (www.eosweb.larc.nasa.gov) for all variables, except UV index data which were from the Tropospheric Emission Monitoring Internet Service (www.temis.nl) as average values of the year 1999. For all variables, we estimated the annual average per geographical locality, and for each species over all localities. We mapped mean values on the species phylogeny and estimated their phylogenetic signal as for morphometrical traits using BayesTraits. We then estimated the coefficient of phylogenetic regression (β) of each morphological trait on each climatic variable using the Continuous model in BayesTraits to determine the sign of the relationship, and estimated its significance by counting the ratio of the time β crossed the zero point. A wellsupported β should minimally switch signs.

In order to be able to obtain an estimate of the "overall effect" of the climatic factors on each trait, we conducted a PCA analysis on the average values of the 14 climatic variables for each species, and then tested the phylogenetic regression using BayesTraits for each quantitative trait and values of the two first principal components, henceforth estimated two scores, denoted ßPC1 and ßPC2, respectively. We then tested the correlation between the absolute values of the regression coefficients, i.e. |ßPC1| and |ßPC2|, with both t and λ for each morphological trait.

Results

Multivariate analysis of morphological traits recovers the two major clades of the melanogaster subgroup

The Bayesian phylogeny of the nine species inferred from the 22 nuclear genes (Figure 1A) recapitulates the current view of relationships based on complete genome sequences of the five annotated species [START_REF] Clark | Evolution of genes and genomes on the Drosophila phylogeny[END_REF], with all nodes having a posterior probability of 100%. The subgroup is divided into two major clades: melanogaster (4 species) and yakuba (5 species).

Principal component analysis (PCA) of the 20 traits in the 11 strains showed a strong concordance with the melanogaster-yakuba phylogenetic split for both females and males (Figure 1B,C). The first four principal components explained 69.5% and 76.1% of the total variance in females and males, respectively. PC1 (34.3% and 34.8% in females and males, respectively) was mostly correlated with abdominal pigmentation traits, PC2 (17.2% and 18.3%) with abdominal bristles in both sexes and ovariole number in females, PC3 (9.8% and 14.5%) with body size, and PC4 (8.2% and 8.5%) with thoracic bristles.

On average (Table S1), the abdomen of the species of the melanogaster clade is three times darker than that of the yakuba clade (summed over the six segments). The two temperate populations of the two cosmopolitan species (D. melanogaster and D. simulans) correctly fall within the morphospace of this clade, indicating the persistence of the pigmentation difference (Figure 1B,C). Within the yakuba clade, D. yakuba is exceptional in having a very bristly abdomen (summed over the six segments) and a dark abdomen comparable to those of the melanogaster clade (Figure 1B,C). The completely yellow abdomen of D. santomea males [START_REF] Lachaise | Evolutionary novelties in islands: Drosophila santomea, a new melanogaster sister species from São Tomé[END_REF] clearly distinguishes it from the remaining species of the subgroup where the last two segments of the male abdomen are completely pigmented (Figure 1C).

Phylogenetic signal of morphological traits does not correlate with short-term heritability

Estimates of the phylogenetic signal (λ) revealed heterogeneity of the effect of the phylogenetic history among trait categories (Table 1). On average over traits and sexes, λ ranged from 0.37 to 0.43 for body size, shape and ovariole number, and from 0.51 to 0.52 for setation and pigmentation traits. There was no correlation between λ in the melanogaster species subgroup and heritability as estimated from the intraclass correlation in D. melanogaster (r = -0.035, P = 0.860) after combining data from Table 1 from both sexes. The lack of correlation also persisted when female (r = -0.122, P = 0.641) and male traits (r = 0.179, P = 0.598) were separately analysed.

Phylogenetic regression reveals a substantial effect of climatic variables on morphological evolution

Principal component analysis (PCA) of the average values of 14 climatic variables deduced from biogeographical data of the nine species of the melanogaster subgroup revealed a strong distinction between the two major clades (Figure 1D). The distinction was most obvious along PC1, which explained 54% of the variance and mostly correlated with solar radiation variables in agreement with the East-West disjunction of the likely origin of the two clades. PC2, however, which explained 35% of the variance, was most correlated with temperature and relative humidity variables. It mostly discriminated among the most closely-related species, i.e. D. simulans from both D. mauritiana and D. sechellia, and D. yakuba from D. santomea.

Phylogenetic signals (λ) for species means of climatic variables fell within two categories:

low (λ = 0.32-0.43) for temperature, air pressure and relative humidity variables, and intermediate to high (λ = 0.53-0.74) for wind speed and radiation variables (Table S2). We found some correspondences between λ estimates for morphological traits and climatic variables and their phylogenetic regressions (ß) (Figure 2). For example, morphological traits with low λ, such as wing length, ovariole number or pigmentation of the last abdominal segment in females (P7) were significantly negatively correlated with climatic variables with low λ such as temperature and relative humidity.

Because ß varied substantially among individual climatic variables, and climatic variables are usually correlated, we conducted the phylogenetic regression analyses on the PC1 and PC2 scores of the species climatic data (Table 1). For all traits and sexes (combined from Table 1), |ßPC1| did not correlate with either t (r = -0.248, P = 0.221) or λ (r = -0.070, P = 0.691). On the contrary, |ßPC2| had a highly significant positive correlation with t (r = 0.585, P = 1.686 x 10 -3 ) and negatively correlated with λ (r = -0.421, P = 0.012).

In summary, short-term heritability estimates do not correlate with phylogenetic signal, i.e. the so-called species-level or phylogenetic heritability, but they do correlate with phylogenetic regression on species-averaged climatic variables related to temperature and aridity.

Discussion

In spite of their conceptual resemblance, i.e. increased similarity due to genetic relatedness, our results do not support a correlation between the measurements of heritability and phylogenetic signal, at least for the studied quantitative morphological traits in the Drosophila melanogaster species group. On the one hand, heritability is explained by the genetic architecture of a trait, whereas phylogenetic signal refers to how evolution in one trait correlates with speciational history.

Low phylogenetic signal indicates high homoplasy, which could be caused by a multitude of evolutionary factors, dependent on or independent from the genetic architecture of the trait. Our observation that quantitative morphological characters in D. melanogaster show a significant phylogenetic regression on climatic variables (Figure 2), suggests that convergent selection on those heritable traits might have been the cause of their low phylogenetic signal, i.e. homoplasy.

Body size, bristle number and pigmentation have been classical traits in Drosophila quantitative genetics for nearly a century, but they have rarely been used in systematics or phylogenetics due to their low phylogenetic signal [START_REF] Sayad | Quantifying the extent of morphological homoplasy: A phylogenetic analysis of 490 characters in Drosophila[END_REF]. In fact, even within D. melanogaster, most of those characters show evidence for convergent evolution among geographical populations inhabiting colder regions (e.g., high altitudes and latitudes on different continents) [START_REF] Gibert | Comparative analysis of morphological traits among Drosophila melanogaster and D. simulans: genetic variability, clines and phenotypic plasticity. In: Drosophila melanogaster, Drosophila simulans: So Similar, So Different[END_REF]). For wing size for example, larger wings in high latitude populations in Australia and South America are mostly driven by different cellular mechanisms, namely cell number proliferation and cell size enlargement, respectively [START_REF] Zwaan | Cellular basis of wing size variation in Drosophila melanogaster : a comparison of latitudinal clines on two continents[END_REF]. In high altitude Ethiopian populations with the largest size traits in D. melanogaster, cell size enlargement was also found to significantly drive the evolution of larger wings (Lack et al. 2016b, c), in spite of the fact that those flies are genetically closer to normal size African populations than to other non-African flies (Lack et al. 2016a). For pigmentation, convergent evolution of melanism occurs between high latitude populations and high altitude African populations, but the underlying genetic basis differs between populations [START_REF] Bastide | Pigmentation in Drosophila melanogaster reaches its maximum in Ethiopia and correlates most strongly with ultra-violet radiation in sub-Saharan Africa[END_REF][START_REF] Bastide | A Variable Genetic Architecture of Melanic Evolution in Drosophila melanogaster[END_REF]. Even when the same genetic locus underlies the same trait within and between species, apparent phylogenetic similarities can still hide distinct genetic origins, such as the light female morph phenotype caused by independent mutations at the same regulatory element in the sister species D. erecta and D. orena [START_REF] Yassin | Ancient balancing selection at tan underlies female colour dimorphism in Drosophila erecta[END_REF].

Our results might have been biased by our choice of characters. Indeed, an analysis of the phylogenetic signal of 490 morphological characters in drosophilids revealed that nearly two thirds of the characters evolve convergently (Al Sayad and Yassin 2019). However, genital traits were among those with the highest signal. In the melanogaster subgroup, male and female genital structures are indeed the most discriminatory characters between species [START_REF] Yassin | Coevolution between male and female genitalia in the Drosophila melanogaster species subgroup[END_REF]. Due to their small size and dissection-associated difficulties, few studies have aimed at investigating their quantitative genetics. It was however estimated that the heritability of the size of the male epandrial posterior lobe, which is used to grasp the female during copulation, was 0.65 in D. melanogaster [START_REF] Takahara | Genome-Wide Association Study on Male Genital Shape and Size in Drosophila melanogaster[END_REF] and as high as 0.93 in crosses between D. simulans and D. mauritiana [START_REF] Zeng | Genetic Architecture of a Morphological Shape Difference Between Two Drosophila Species[END_REF]. Here, a correspondence between high heritability and phylogenetic signal for the same trait may be expected, if the trait evolution has contributed to the phylogenetic divergence between species.

In conclusion, our results caution against the usage of the term heritability for phylogenetic signal, even if some precisions such as phylogenetic or species-level heritability are applied.

Phylogenetic signal is a description of a pattern that can be generated by a multitude of processes [START_REF] Losos | Seeing the Forest for the Trees: The Limitations of Phylogenies in Comparative Biology[END_REF]. Although a correspondence between the two concepts may sometimes be found, mostly at short time scales (e.g., in viruses, [START_REF] Shirreff | How effectively can HIV phylogenies be used to measure heritability?[END_REF][START_REF] Vrancken | Simultaneously estimating evolutionary history and repeated traits phylogenetic signal: applications to viral and host phenotypic evolution[END_REF], it is most likely, given the preponderance of homoplasy for most types of characters, that a strict correlation will be rare. Our results also demonstrate a potential role of environmental conditions in shaping the genetic architecture of a trait. The positive correlation between short-term heritability and long-term phylogenetic regression on variables related to temperature and aridity may be due to the fact intraclass correlations estimates of heritability may either involve a part of the genotype x environmental variance (i.e. broad-sense heritability) or to the direct action of selection on the genetic variance underlying quantitative traits. Future dissection of the genetic architecture of the traits analysed here in the remaining species of the melanogaster subgroup would definitively help explaining the role of the environment in the long-term evolution of those classical Drosophila quantitative traits.

  simulans, we used two strains, one from a temperate origin, Prunay (France), collected in 2010 and the other from the island of Sao Tomé at a latitude close to the equator, provided by D. Lachaise in 2002. D. sechellia was founded in 1985 by mixing 25 isofemale lines. D. mauritiana was founded in the same way in 1988. For D. teissieri, we used the type strain, established in 1970 from Mount Selinda in Zimbabwe. For D. yakuba we used a still older strain, collected in Kunden (Cameroon) in 1967. For D. santomea, we used a mass culture given by D. Lachaise and established in 2005 from Sao Tomé. The last two species corresponded to isofemale lines. For D. erecta, we used a line collected in 2005 by D. Lachaise in the park of La Lopé (Gabon), and homozygous for the female light abdomen locus. D. orena was the only line available so far, collected in the mountains of Cameroon in 1975.

Figure 1 -

 1 Figure 1 -(A) Bayesian phylogeny of the nine species of the melanogaster subgroup inferred from

Figure 2 -

 2 Figure 2 -Heatmap of (1 -P) where P is the significance probability of the phylogenetic regression

Table 1 -

 1 Phylogenetic signal (λ), heritability as intra-class correlation among isofemale lines (t) and average absolute coefficient of phylogenetic regression |ß| of PC1 and PC2 scores of climatic variables for 44 morphological traits in the nine species of the Drosophila melanogaster species subgroup. Trait abbreviations are given in Methods.

	Trait	Females	Males
		t	λ	|ßPC1| |ßPC2| t	λ	|ßPC1| |ßPC2|
	Body size
	WL	0.50 0.40 0.82 0.98 0.44 0.44 0.77 0.96
	TL	0.38 0.33 0.92 0.56 0.29 0.38 0.96 0.68
	TW		0.43 0.92 0.91	0.46 0.94 0.96
	Body shape
	WL/TL	0.37 0.46	0.37 0.44
	WL/TW		0.48		0.47
	TL/TW		0.36		0.37
	Setation
	STPn	0.31 0.73 0.92 0.52 0.28 0.53 0.77 0.92
	STPa		0.76 0.99 0.75	0.61 0.96 0.92
	STPp		0.37 0.91 0.98	0.38 0.94 0.96
	STPr		0.79		0.62
	A2	0.11 0.31 0.74 0.73 0.21 0.32 0.82 0.71
	A3	0.19 0.60 0.79 0.52 0.23 0.42 0.64 0.89
	A4	0.23 0.63 0.68 0.61 0.15 0.45 0.72 0.68
	A5	0.12 0.64 0.63 0.51 0.06 0.40 0.75 0.63
	A6	0.09 0.36 0.59 0.57
	A7	0.17 0.37 0.76 0.68
	T7		0.46 0.97 0.89
	Pigmentation
	P2	0.26 0.65 0.75 0.57 0.26 0.42 0.58 0.74
	P3	0.21 0.69 0.70 0.57 0.28 0.51 0.50 0.64
	P4	0.22 0.68 0.60 0.53 0.23 0.54 0.65 0.69
	P5	0.36 0.62 0.71 0.71	0.40 0.75 0.81
	P6	0.41 0.49 0.55 0.92	0.44 0.77 0.83
	P7	0.33 0.41 0.65 0.99
	Reproductive organs
	OV	0.37 0.37 0.88 0.99

Conflicts of interest

We declare no conflict of interest.


Angela Villota

Diaz Jesus Aranda

Juan Francisco

Juan Francisco Modelando

Ángela Villota

Juan Francisco Díaz

Modelando Sistemas de Membranas en ntcc

Keywords: ntcc, membrane systems, systems biology

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introducción

Los sistemas biológicos pueden ser vistos como sistemas en donde existen procesos que interactúan de forma concurrente y evolucionan en el tiempo. El área que aborda la biología desde esta perspectiva es llamada biología de sistemas [START_REF] Kitano | Systems Biology: A Brief Overview[END_REF][START_REF] Kitano | Systems Biology: Toward System-level Understanding of Biological Systems[END_REF]. En el contexto de la biología de sistemas se ha evidenciado la utilidad de las herramientas computacionales que permiten modelar y simular procesos biológicos conocidos. Entre estas herramientas se encuentran los cálculos de procesos que han mostrado ser útiles para modelar, simular y verificar propiedades de sistemas biológicos [START_REF] Gutierrez | Modelamiento de Sistemas Biologicos usando Calculos de Procesos Concurrentes[END_REF][START_REF] Gutierrez | Timed Concurrent Constraint Programming for Analysing Biological Systems[END_REF][START_REF] Busi | Modeling and Analysis of Biological Processes by Mem(Brane) Calculi and Systems[END_REF][START_REF] Regev | BioAmbients: an abstraction for biological compartments[END_REF][START_REF] Cardelli | Brane Calculi Interactions of Biological Membranes[END_REF][START_REF] Arbelaez | A Generic Framework to Model, Simulate and Verify Genetic Regulatory Networks (Poster)[END_REF]. De entre las distintas posibilidades, este trabajo explora la expresividad de ntcc, un cálculo basado en el modelo CCP [START_REF] Saraswat | The Semantic Foundatiosn of Concurrent constraint Programming[END_REF] (Concurrent Constraint Programming), en particular su adecuación para modelar sistemas de membranas. En [START_REF] Gutierrez | Modelamiento de Sistemas Biologicos usando Calculos de Procesos Concurrentes[END_REF][START_REF] Gutierrez | Timed Concurrent Constraint Programming for Analysing Biological Systems[END_REF][START_REF] Arbelaez | A Generic Framework to Model, Simulate and Verify Genetic Regulatory Networks (Poster)[END_REF] se muestra la utilización de ntcc para modelar procesos biológicos destacando los beneficios que trae la representación de información parcial, la independencia entre los constructores del cálculo y la demostración de propiedades a través de una lógica temporal (LTL).

Por otra parte, existen formalismos (Brane [START_REF] Cardelli | Brane Calculi Interactions of Biological Membranes[END_REF], BioAmbients [START_REF] Busi | Modeling and Analysis of Biological Processes by Mem(Brane) Calculi and Systems[END_REF], Membrane Systems [START_REF] Paun | Computing with Membranes[END_REF]) que modelan el comportamiento de membranas biológicas y que son una alternativa diferente y sencilla para la descripción y simulación de sistemas biológicos. En este artículo se propone un método que permite describir sistemas biológicos en ntcc de forma similar a los formalismos que modelan membranas (Sistemas de Membranas). Inicialmente se realizó una comparación de la forma en que se modelan sistemas en Brane y BioAmbients, con respecto a ntcc, y en general con cualquier derivado de CCP. Se encontraron tres diferencias importantes:

El mecanismo de comunicación entre procesos: en CCP la comunicación es asíncrona y se hace por medio del store que es visible para todos los procesos. Los estímulos o entradas están dados por medio de restricciones que los procesos pueden imponer o consultar. En BioAmbients/Brane, la comunicación es síncrona y requiere de la presencia de acciones y coacciones, o de canales de comunicación en los procesos. La forma en que se construyen los sistemas: en BioAmbients/Brane las membranas son contenedores de procesos que pueden tener otras membranas al interior (anidamiento), por tanto los sistemas en estos cálculos pueden ser vistos como jerarquías de membranas. En CCP los sistemas pueden ser especificados de forma modular, pero no existen jerarquías entre ellos. Movilidad: los cálculos de membranas proveen un conjunto de operaciones que permiten modelar el ingreso/salida de membranas; esto puede realizarse gracias a que existe movilidad como concepto del cálculo. En CCP no hay movilidad como concepto del cálculo y aunque es posible modelar la movilidad como está propuesto en [START_REF] Laneve | Mobility in the cc-paradigm[END_REF], no es posible modelar el mismo comportamiento que en BioAmbients y Brane.

Las anteriores diferencias explican las dificultades que se presentan para modelar las operaciones de membranas con los constructores de ntcc. Pero esto no quiere decir que no se puedan modelar sistemas biológicos teniendo en cuenta la localización de sustancias y procesos (al estilo de membranas). Para ello es necesario definir un conjunto de procesos y variables que hagan posible modelar el comportamiento de las operaciones de membranas y especificar sistemas en los cuales exista una jerarquía de membranas en ntcc.

Este documento está organizado así: en la sección 2 se encuentra una breve descripción de la sintaxis de ntcc. En la sección 3 se describen las técnicas que permiten modelar sistemas de membranas en ntcc. En la sección 4 se aplican las técnicas de la sección 3 para modelar la ruta endocítica de la LDL del colesterol. La sección 5 presenta las conclusiones y el trabajo futuro.

El Cálculo ntcc

El cálculo de programación concurrente por restricciones temporal, ntcc [START_REF] Valencia | Temporal Concurrent Constraint Programming[END_REF] es un formalismo para sistemas reactivos. La noción de tiempo es introducida por medio de una secuencia de unidades de tiempo, en donde cada unidad tiene su propio almacén de restricciones (store) y dura mientras los procesos realizan sus computaciones y añaden (tell ) o infieren (ask ) información del store.

Los procesos ntcc están parametrizados en un sistema de restricciones, el cual especifica qué tipo de restricciones puede manejar el store. Los procesos P, Q, . . . ∈ P roc son construidos como se muestra en el cuadro 1.

El proceso tell(c) adiciona la restricción c al store, dejando a c disponible para los otros procesos en la unidad de tiempo actual. La selección con guarda i∈I when c i do P i , donde I es un conjunto finito de índices, representa un proceso que, en la unidad de tiempo actual, escoge de manera no determinista un proceso P j (j ∈ I) cuya guarda c j es deducida del store (las alternativas no escogidas son descartadas). El proceso P Q representa la composición paralela entre P y Q (en una unidad de tiempo P y Q operan concurrentemente, comunicándose por medio del store). El proceso (local x) P se comporta como P pero la información de la variable x es local a P , i.e. P no puede ver la información sobre una variable global x y los procesos que no hacen parte de P no pueden ver la información producida por P sobre x. El proceso next P ejecuta el proceso P en la siguiente unidad de tiempo. El proceso unless c next P representa la activación de P en la siguiente unidad de tiempo si y solo si c no puede ser deducida a partir del store en la unidad de tiempo actual. El proceso !P ejecuta P en todas las unidades de tiempo a partir de la actual (esta es la versión de retraso del operador replicación del cálculo π). El proceso ⋆P representa un retraso sin límite pero finito de P , i.e. P finalmente será ejecutado en alguna unidad de tiempo en el futuro.

P, Q, . . . ::= tell(c) Tell | P i∈I when ci do Pi Selección con guarda | P Q Composición Paralela | (local x) P Comportamiento Local | next P Retraso de Unidad | unless c next P Vencimiento | !P Replicación | ⋆P Retraso acotado
Cuadro 1. Sintaxis de ntcc

Técnicas para la Modelación de Sistemas de Membranas en ntcc

A continuación se presentan dos técnicas que son claves para modelar sistemas de membranas: modelar membranas, su localización y jerarquía (ver sección 3.1), y modelar reacciones químicas, las cuales son el mecanismo de interacción entre los objetos situados en las membranas (ver sección 3.2).

Modelando Membranas, Localización y Jerarquías

Los procesos y definiciones propuestos a continuación están inspirados en la forma en que se modelan las membranas con Sistemas de Membranas (P Systems) que aparecen en [START_REF] Busi | Modeling and Analysis of Biological Processes by Mem(Brane) Calculi and Systems[END_REF].

Para la especificación de sistemas al estilo de membranas se utilizará un sistema de restricciones sobre dominios finitos. La ubicación de los elementos del sistema se realizará por medio de variables. La estructura del sistema es estática, los contenidos de las membranas son modificados debido al paso de sustancias a través de las membranas (entrada y salida de sustancias), pero no hay desplazamiento de membranas como tal. Una membrana en este enfoque, es una frontera que encierra una región llamada nivel. Los niveles que conforman un modelo pueden variar dependiendo de la descripción del proceso biológico que se desea modelar.

Al interior de cada nivel se encuentra ubicado un conjunto de objetos que representan las sustancias presentes en la descripción biológica del sistema. La ubicación de cada uno de estos objetos (y de los compuestos que conforman), está dada por una variable declarativa que puede tomar valores en los naturales, y denota la cantidad de cierto objeto (y compuesto) en un nivel dado.

En el cuadro 2 se muestra un ejemplo sencillo de un sistema con dos niveles, en donde se han ubicado tres tipos de objetos A, B y C. Por cada tipo de objeto (v.gr. A) se crean dos variables (A 0 y A 1 , en el caso de A). Los valores de dichas variables representan la presencia y cantidad de los objetos de un tipo dado en cada nivel.

nivel 0 C nivel 1 A, A, B Variables del nivel 0: A0 = 0, B0 = 0, C0 = 1 Variables del nivel 1: A1 = 2, B1 = 1, C1 = 0
Cuadro 2. Ejemplo: Sistema con dos niveles y cuatro objetos.

Modelando Reacciones Químicas

Las reacciones químicas son procesos en los cuales una o más sustancias (reactivos) sufren transformaciones químicas para convertirse en otras (productos). Las reacciones químicas que se modelan a continuación son reacciones de síntesis, descomposición y transporte, o combinaciones de ellas. Una reacción de síntesis es aquella en que dos o más sustancias se unen para formar un compuesto. En una descomposición, una sustancia compleja (varias sustancias enlazadas) se divide en sus componentes. Las reacciones de transporte permiten mover una sustancia a través de las membranas del sistema.

A, B → AB (síntesis), AB → A, B (descomposición), [A]1 → [ ]1, A (transporte)
De forma general, una reacción química de la forma R1, . . . , Rn → P1, . . . , Pm es modelada en ntcc con el proceso Reaction. El proceso Reaction es genérico y permite modelar los tipos de reacciones mencionados (ver cuadro 3). Con el fin de dar continuidad a la evolución del proceso, los valores de las variables que están por fuera de una reacción son pasados a la siguiente unidad de tiempo sin alteración; esto es representado por Rest. Para evitar inconsistencias es necesario que solo una reacción pueda ser activada en cada unidad de tiempo; esto puede lograrse usando selección no determinista (+) entre las reacciones y niveles.

Reaction def = when (R1 ≥ 1 . . . ∧ Rn ≥ 1) do next ( tell( R1 = R1 -1 ∧ . . . ∧Rn = Rn -1 P1 = P1 + 1 ∧ . . . ∧ Pm = Pm + 1 ∧ Rest) )
Nótese que el antecedente de una reacción es un conjunto de restricciones que pueden ser deducidas (o no) del store. En el caso en que estas restricciones puedan ser deducidas, se imponen otras en la siguiente unidad de tiempo (esto sucede a causa del constructor next).

En adelante los procesos que modelen reacciones químicas y que sean derivados de Reaction se representarán así: condiciones → consecuencia.

when (A ≥ 1 ∧ B ≥ 1) do next (A&B = A&B + 1 ∧ A = A -1 ∧ B = B -1) (1) when (A&B ≥ 1) do next (A = A + 1 ∧ B = B + 1 ∧ A&B = A&B -1) (2)
when En la primera parte de esta sección se encuentra una descripción del sistema biológico que se desea modelar. En la segunda parte se presenta el modelo propuesto usando las técnicas presentadas en la sección anterior.

(Ai ≥ 1) do next (Aj = Aj + 1 ∧ Ai = Ai -1) (3)

Ruta Endocítica de la LDL del Colesterol

El conjunto de pasos por medio de los cuales una célula lleva a su interior una partícula, para luego ser procesada por sus organelas, es conocido como ruta endocítica. Las células llevan a su interior partículas de colesterol siguiendo una ruta endocítica. El colesterol puede ser de dos clases: de alta densidad (high density) o de baja densidad (low density). El colesterol de baja densidad se transporta por el torrente sanguíneo al interior de una lipo-proteína conocida como la LDL (low density lipoprotein). En la Figura 1 (tomada de [START_REF] Lodish | Molecular cell biology[END_REF]) se muestra el proceso de la ruta endocítica de la LDL del colesterol.

El primero de los pasos en la ruta endocítica se llama endocitosis mediada por receptor. El proceso se inicia cuando un receptor (ubicado en la membrana celular) enlaza a la partícula de LDL (ligando), a través de la proteína apoB100 que está ubicada en la superficie de la partícula (paso 1 de la figura 1).

Una vez ocurrido el enlace, el complejo receptor-partícula ingresa a la célula en una vesícula de transporte formada por endocitosis; ésta vesícula está rodeada por clatrina (proteína que recubre la membrana). En el momento en que la clatrina se desprende de la vesícula, ésta se convierte en un endosoma temprano (pasos 2 y 3 de la figura 1).

A continuación, el endosoma temprano se fusiona con un endosoma tardío (una de las organelas de la célula). El ambiente al interior del endosoma tardío es ácido en comparación con el exterior (ph entre 5 y 6), lo que causa que se rompa el enlace entre el receptor y el ligando. Después de la separación, el receptor debe ser devuelto a la superficie de la membrana celular con el fin de servir como mecanismo para la captura de otras partículas de LDL al exterior de la célula. El receptor es enviado en el interior de una vesícula que se desprende del endosoma tardío (paso 5 de la figura 1).

El endosoma tardío que contiene únicamente la partícula de LDL se fusiona con un lisosoma. El lisosoma es una organela que contiene las enzimas necesarias para degradar la cubierta de la partícula y extraer el colesterol de su interior para luego ser aprovechado por la célula (paso 6 de la figura 1).

Figura 1. Ruta endocítica de la lipoproteína de baja densidad del colesterol, tomada de [START_REF] Lodish | Molecular cell biology[END_REF]

Modelo en ntcc

Con base en la descripción biológica de [START_REF] Lodish | Molecular cell biology[END_REF] y descrita anteriormente, se proponen cuatro niveles como estructura del sistema. El conjunto de objetos que van a localizarse en los niveles definidos anteriormente es el siguiente: S = {LDLpart, LDLRcpt, Clth, LateEndo} en donde LDLpart, representa la lipoproteína, LDLRcpt representa el receptor, Clth representa a la proteína que cubre la vesícula que ingresa a la célula (clatrina) y LateEndo es el endosoma tardío. Los objetos anteriormente mencionados pueden formar compuestos. Estos compuestos se representan por la unión de los nombres de los objetos con un "&" intermedio; por ejemplo LDLpart&LDLRcpt representa al compuesto que resulta al unirse un ligando con el receptor. Cada uno de los niveles modela una región de la célula en donde hay sustancias e interacción escrita como procesos ntcc, que agrupan las reacciones que pueden ocurrir en cada uno de los niveles.

Nivel 0. Modela el exterior de la célula, en este nivel está localizada la lipoproteína de LDL. No hay reacciones en este nivel.

Nivel 1. Modela la membrana celular, no como una división, sino como el espacio en donde está localizado el receptor. En esta región se localiza una reacción que permite la unión entre el ligando (LDLpart) y el receptor (LDLRcpt), así como también la localización del compuesto formado en el nivel 1 (membrana de la célula).

Nivel 2. Modela el interior de la célula, aquí están localizados el endosoma tardío y el lisosoma. Comprende seis reacciones que modelan el comportamiento al interior de la célula. El proceso R 2,1 modela el paso 2 de la figura 1, en donde una vesícula envuelta en clatrina que contiene el receptor y el ligando unidos ingresa a la célula. Por medio de los procesos R 2,2 y R 2,3 se modela el momento en que la clatrina se desprende de la vesícula, convirtiéndose en endosoma temprano y la fusión de éste con el endosoma tardío. En R 2,4 se modela la ruptura del enlace entre el receptor y el ligando, posteriormente en R 2,5 se modela el paso del receptor del nivel 2, al nivel 1 para que sea reutilizado, como se puede ver en el paso 5 de la figura 1. Por medio del proceso R 2,6 se modela la fusión entre el endosoma tardío y el lisosoma. Nivel 3. Modela el interior del lisosoma; este es el nivel a donde debe llegar la partícula de LDL para ser degradada por la célula. Este nivel alberga tres reacciones, los procesos R 3,1 y R 3,3 modelan el procesamiento de la LDL del colesterol y el proceso R 3,2 permite retornar el objeto LateEndo al nivel 2, para que pueda ser reutilizado en un nuevo proceso.

Procesos adicionales. Los anteriores niveles componen el proceso CellActivity por medio de + (selección no determinista) para evitar que dos o más reacciones de distintos niveles, sean activadas y causen inconsistencias. La característica principal de este proceso es la persistencia. Por medio del operador ! se garantiza que durante la evolución del sistema, los procesos que simulan la célula están disponibles en todas las unidades de tiempo.

Init es un proceso que asigna valores iniciales a las variables otorgando una configuración inicial para la simulación; en este caso en particular se ubican un lisosoma, y un endosoma temprano en el nivel 2 (interior de la célula) y un receptor sobre la membrana (nivel 1). El proceso principal es System, el cual contiene la actividad de la célula en paralelo con las condiciones iniciales y el proceso tell(LDLpart 0 = n 0 ) que inicia la evolución del proceso, con n 0 partículas de LDL del colesterol al exterior de la célula. Todos los procesos aquí descritos se encuentran en el cuadro 4.

Level1 def = LDLpart0, LDLRcpt1 → LDLpart&LDLRcpt1 Level2 def = R2,1 + R2,2 + R2,3 + R2,4 + R2,5 + R2,6 R2,1 = Clth2, LDLpart&LDLRcpt1 → Clth&LDLpart&LDLRcpt2 R2,2 = Clth&LDLpart&LDLRcpt2 → Clth2, LDLpart&LDLRcpt2 R2,3 = LateEndo2, LDLpart&LDLRcpt2 → LateEndo&LDLpart&LDLRcpt2 R2,4 = LateEndo&LDLpart&LDLRcpt2 → LateEndo&LDLpart2, LDLrcpt2 R2,5 = LDLrcpt2 → LDLrcpt1 R2,6 = LateEndo&LDLpart2 → LateEndo&LDLpart3 Level3 def = R3,1 + R3,1 + R3,1 R3,1 = LateEndo&LDLpart3 → LateEndo3, LDLpart3 R3,2 = LateEndo3 → LateEndo2 R3,3 = LDLpart3 → Chol3 CellActivity def = !( Level 1 + Level 2 + Level 3 ) Init def = tell (LDLrcpt1 = 1 ∧ Clth2 = 1 ∧ LateEndo2 = 1) System def = tell (LDLpart0 = n0) Init CellActivity
Cuadro 4. Modelo de la ruta endocítica de la LDL del colésterol con el cálculo ntcc La selección no determinista entre las reacciones y los niveles es útil porque permite garantizar que solo una reacción pueda ser activada en cada unidad de tiempo; esto es necesario ya que las variables declarativas no cambian su valor. No obstante el modelo no garantiza la ejecución de cada proceso (activación de una reacción) en un momento dado, e.g con el proceso !(r 1 + r 2 ) no se puede garantizar que en alguna unidad de tiempo r 1 sea ejecutado. A continuación se muestra un segundo modelo que incluye el operador entre las reacciones en cada nivel y los niveles en CellActivity. La composición paralela de reacciones causa inconsistencias, porque permite la activación de varias de ellas en la misma unidad de tiempo; por tanto también es necesario modificar el proceso Reaction e incluir el proceso Control i que evita que los valores de las variables que hacen parte de la reacción sean modificados por otros procesos. En CellActivity se incluye el proceso Control gen que se encarga de actualizar los valores de las variables que no hacen parte de una reacción y preservar la continuidad de la evolución del sistema. El proceso Control gen se compone de varios procesos en paralelo, uno por cada variable. Dichos procesos son de la forma:

unless Control_var = 1 do next var ∀ var ∈ V ars, donde V ars es el conjunto de variables del sistema. Para el funcionamiento de los procesos de control se declara una variable de control por cada una de las variables del sistema. El nombre de cada variable de control está dado por el prefijo Control_ y el nombre de la variable con la que está relacionada. Por ejemplo para la variable LDLpart 0 la variable de control correspondiente tiene el nombre Control_LDLpart 0 . En el cuadro 5 se encuentran los procesos que cambian para el segundo modelo.

Reaction def = when (R1 ≥ 1 . . . ∧ Rn ≥ 1) do next ( tell( R1 = R1 -1 ∧ . . . ∧Rn = Rn -1 P1 = P1 + 1 ∧ . . . ∧ Pm = Pm + 1) ) Controli Controli def = Control_Ri = 1, Control_P j = 1 para 1 ≤ i ≤ n y 1 ≤ j ≤ m Level1 def = LDLpart0, LDLRcpt1 → LDLpart&LDLRcpt1 Level2 def = R2,1 R2,2 R2, 3

Conclusiones y Trabajo Futuro

Con este trabajo se muestra que es posible modelar sistemas de membranas en ntcc. Las técnicas propuestas pueden ser aplicadas para modelar procesos biológicos que consideren la interacción de membranas. Sin embargo, existen aún inconvenientes con la selección no determinista y la composición paralela entre las reacciones y los niveles.

El problema que causa la selección no determinista interfiere con la automatización de la simulación del proceso biológico, ya que no se puede garantizar que la máquina alcance algunos estados del proceso. Este problema es resuelto con la composición paralela, pero requiere de un cuidadoso manejo de las variables declarativas. Una variable declarativa no puede cambiar su valor, por tanto no se puede permitir que dos reacciones que involucran una misma sustancia sean activadas en la misma unidad de tiempo. Por ejemplo, en el caso en que existan más de una partícula de LDL al exterior de la célula y más de un receptor en la membrana, se genera una inconsistencia cuando se imponen las restricciones LDL_part&LDL_Rcpt = LDL_part&LDL_Rcpt + 1, y LDL_part&LDL_Rcpt = LDL_part&LDL_Rcpt -1 en la misma unidad de tiempo. La inconsistencia no se presenta cuando hay más de una partícula de LDL al exterior y un solo receptor, así que en este caso la competencia por un recurso no es un problema. La utilización de composición paralela hace que el modelo sea más cercano al proceso biológico.

Las técnicas propuestas en este trabajo modelan membranas como fronteras de los niveles. También es posible modelar interacción al interior de una membrana, como ocurre en la ruta endocítica de la LDL del colesterol. Para lograr esto, se incluyó en el modelo un nivel (nivel 1) y se ubicó en dicho nivel los objetos y reacciones que representan el comportamiento de la membrana celular en dicho proceso biológico.

Todo apunta a que la mejor alternativa es conservar la composición paralela entre las reacciones y los niveles, pues así el modelo es más fiel al sistema biológico, y es posible automatizar la simulación y verificar propiedades. Las inconsistencias son un problema inherente al tipo de variable y los dominios. Una propuesta interesante es examinar la utilización de otros sistemas de restricciones, como por ejemplo flujos, en los cuales se puedan actualizar los valores de las variables; así desaparecerían las inconsistencias. Como trabajo futuro se propone explorar sistemas de restricciones distintos a los dominios finitos, con el fin de dar solución al problema de las inconsistencias. Adicionalmente, se espera aplicar la lógica temporal (LTL) para verificar propiedades en el modelo.

Cuadro 3 . 4 .

 34 Reacciones químicas en ntcc. (1): Síntesis, las sustancias A, B se unen para formar el compuesto A&B. (2): Descomposición, el compuesto A&B se separa en las sustancias A, B. (3): Transporte, la sustancia A se mueve de un nivel i a un nivel j Modelo de la Ruta Endocítica de la LDL del Colesterol en el Cálculo ntcc

Level 1

 1 Level 2 Level 3 Controlgen) Controlgen def = Π unless Control_var = 1 do next var, ∀ var ∈ V ars Cuadro 5. Modelo con paralelismo y el proceso Controlgen

* Research partially supported by the INRIA Equipe Associee FORCES. The work of Jesús Aranda has been supported by COLCIENCIAS (Instituto Colombiano para el Desarrollo de la Ciencia y la Tecnologia "Francisco Jose de Caldas"), INRIA and EGIDE (Centre francais pour l' accueil et les echanges internationaux).


Owen Gauthier
email: owen.gauthier@grenoble-inp.org

Sébastien Haendler

Ronan Beucher

Patrick Scheer

Quentin Rafhay

Christoforos Theodorou

Static and LFN/RTN Local and Global Variability Analysis Using an Addressable Array Test Structure

Keywords: matching, low frequency noise, Random Telegraph Noise, variability, addressable test structure, statistical analysis

The use of an addressable array test structure designed on a 28 nm FD-SOI technology for the variability analysis of static, low frequency noise (LFN) and Random Telegraph Noise (RTN) matching is presented. The experimental setup was validated, and a statistical analysis of the above electrical quantities is provided. Using such structures, combined with a switching matrix, local and global variability analysis can be performed while significantly increasing the number of samples, thus enabling a better description of the variations in LFN and RTN, especially when RTN signatures can be scarce. We show that local variations dominate the noise variability compared to global variations.

I. INTRODUCTION

The reduction of the device area in CMOS technologies leads to greater variability of transistor matching, low frequency noise (LFN) and Random Telegraph Noise (RTN), which seriously affects the performances of analog and digital circuits [START_REF] Takeuchi | Single-charge-based modeling of transistor characteristics fluctuations based on statistical measurement of RTN amplitude[END_REF]. Due to its intrinsic stochastic nature, a significant number of samples is needed to properly assess the LFN and RTN parameters variability, especially for small area devices, knowing that LFN and RTN variability scales with the inverse of the device area [START_REF] Ghibaudo | Impact of Scaling Down on Low Frequency Noise in Silicon MOS Transistors[END_REF], [START_REF] Hung | Random telegraph noise of deep-submicrometer MOSFET's[END_REF]. To this end, the use of an addressable array test structure is proposed in this work to allow the characterization of multiple devices of different sizes from a 28 nm FD-SOI technology [START_REF] Planes | 28nm FDSOI technology platform for high-speed low-voltage digital applications[END_REF]. The test structure design enables the study of "global" variability corresponding to variations across the wafer, and "local" variability corresponding to variations on a single die at the minimal distances available. The test structure is first described, followed by a statistical analysis of static and noise parameters which validates the measurements and offers new results on the variability of these quantities on local and global scales.

II. TEST STRUCTURE DESCRIPTION

A. Test Structure General Design

The addressable test structure designed for this work, represented in Fig. 1, consists of an array of 2048 DUTs, arranged in 2 columns of 1024 identically designed and close pairs, that allows the measuring of the threshold voltage mismatch σ∆V t , similarly to [START_REF] Cros | Analysis of process impact on local variability thanks to addressable transistors arrays[END_REF]. The test structure is embedded as a usual parametric test scribe using 25 pads with a 80 µm pitch,

Matched pair

Geometry with the entirety of it fitting in a 0.4 mm 2 region under the test scribe. A power supply is needed for the operation of the decoder, where the positive power supply voltage is set at V DD = 1.8 V and the negative power supply V SS is set at ground.

B. Test Structure Operation

Each DUT terminal is selected by switches powered by the internal circuit as shown in Fig. 2. The switch transistors consist of pairs of NMOS and PMOS of width W = 4 µm, length L = 0.15 µm and W = 12 µm, L = 0.15 µm respectively. They are designed and placed to operate at a very low voltage and have an on-state resistance R ON ≈ 130 Ω, which will limit the maximum current where the voltage drop across the switch remains negligible. Kelvin accesses are available on the drain, gate, and source terminals for each DUT, and the drain of the unselected DUTs is set with a dedicated pin at the same value as the drain voltage of the DUT to further minimize the leakage of drain access switches.

A 11-bit decoder, part of the addressing circuitry, addresses the 2048 transistors, where the first bit of address corresponds to the pair's left DUT if "0" or right DUT if "1". The following bits select the row of the transistor array.

C. Experimental Setup and Measured Devices

The experimental setup is divided in two parts: on one side the noise measurement setup, and on the other side the V DD is set as input. Thus, the 11 bits can be efficiently driven with one single voltage power supply. All the nodes (DUT terminals and addresses) are then connected to a 25-pin probe card mounted on an industrial prober.

V g V bulk V g V s V s V d off V d V
To evaluate the viability of the test structure, NMOS transistors of EOT = 14 Å and a low V t flavor were tested. A mismatch analysis was performed on 16 different geometries. For the noise analysis, the minimal accessible area transistor of W = 0.08 µm, L = 0.03 µm and a larger one of W = 0.5 µm, L = 0.5 µm were chosen. DC drain current (I d ) and drain current noise (S I d ) power spectral density (PSD) measurements were performed for various gate voltages V g bias in the linear regime (V d = 0.05 V). The studies include global variability analysis corresponding to the statistical study across the wafer with one measured device per die for around 60 dies, and local variability analysis on a single die level with 100 measured devices on one die. Measurements of those devices on standard test structures, independent of the one designed for this work, were also performed as references to validate the results obtained from the addressable test structure. device, respectively. The current values and variability are a first indicator of the viability of the test structure, indicating that the measured currents correspond to the addressed devices and not to surrounding switch transistors. This is further confirmed by the measurements of devices of the same shapes from the independent standard test structure which exhibit the same median characteristics as the ones measured with the addressable test structure, as shown Fig. 3. Finally, the drain voltage across the switch transistors was measured by taking advantage of the Kelvin accesses of the DUT, showing a voltage drop of less than 5 mV under a gate voltage operation of 1 V, and thus demonstrating their negligible impact on the measurement of electrical characteristics of the DUT.

III. RESULTS AND DISCUSSION

A. Static Mismatch Analysis

The single die mismatch analysis of the threshold voltage, extracted with the constant current method, is shown Fig. 4, where the extracted values clearly follow Pelgrom's law [START_REF] Pelgrom | Matching properties of MOS transistors[END_REF] with a coefficient A ∆Vt = 1.3 mV.µm. The test structure allowed a comparison to be conducted between an analysis within a single die for 64 matching pairs per geometry, and between all 60 dies of a wafer, showing similar results and dispersion Fig. 5, demonstrating that local variations are dominant in the static variability processes.

B. LFN Variability Analysis

The LFN variability analysis is undertaken using the inputreferred equivalent gate voltage noise PSD S Vg = S I d /g 2 m , g m being the transconductance, because it is suppressing the DC-related variability assuming that LFN is dominated by charge trapping/detrapping [START_REF] Theodorou | Low frequency noise variability in ultra scaled FD-SOI n-MOSFETs: Dependence on gate bias, frequency and temperature[END_REF]. Examples of the noise spectra global spread of device #1 of the large and smaller geometries at a selected gate voltage V g = 0.4 V are shown in Fig. 6, including their median and extremum values. The median normalized drain current noise S I d /I 2 d extracted at f = 10 Hz from the global samples for the two geometries are shown Fig. 7 with the carrier number fluctuations with correlated mobility fluctuations (CNF/CMF) model fit, given by [START_REF] Hung | A unified model for the flicker noise in metal-oxide-semiconductor field-effect transistors[END_REF], [START_REF] Ghibaudo | Improved Analysis of Low Frequency Noise in Field-Effect MOS Transistors[END_REF]:

S I d I 2 d = g m I d 2 S V f b 1 + Ω I d g m 2 (1)
In (1), S V f b is the flat-band voltage spectral density given by:

S V f b = q 2 kT N t λ t W L C 2 ox f (2)
where N t is the oxide trap density, λ t is the tunnel attenuation distance, and Ω = α sc µ eff C ox , with α sc the Coulomb scattering coefficient. The CNF/CMF fit allows the extraction of the oxide trap density, with N t = 3.0 • 10 17 cm -3 .eV -1 for the large geometry and N t = 2.4 • 10 17 cm -3 .eV -1 for the small geometry. The global median S Vg extracted at f = 10 Hz for the two geometries is then compared in Fig. 8 to the local median S Vg for 100 devices of one die and the reference values obtained from the independent standard structure, showing identical noise levels and thus the same extracted oxide trap density. The variability analysis is carried out by considering the standard deviation of log(S Vg ) since the measured data follows a log-normal distribution [START_REF] Ioannidis | Impact of low-frequency noise variability on statistical parameter extraction in ultra-scaled CMOS devices[END_REF]. Fig. 9 shows global and local σ[log(S Vg )] extracted at f = 10 Hz versus the gate voltage for the two different geometries. A slight increase in standard deviation is observed when increasing the gate voltage, which could be related to the CMF impact on the noise in strong inversion as seen in Fig. 8. The standard deviation of log(S Vg ) was also measured for additional geometries to allow the study versus the inverse square root of device area Fig. 10. A 0.2 decade difference in standard deviation between the global and local samples of the small geometry is observed but can be explained by the stochastic nature of the choice of the samples, leading to the same conclusion as for the static mismatch variability: local variations seem to dominate the noise variability processes. With the same noise levels and variability measured whether the sample is global with devices across the wafer, or local with multiple devices on a single die, the use of the addressable structure to perform statistical analysis of LFN is validated. RTN analysis with numerous devices can then be undertaken.

S

C. RTN Statistical Analysis

A statistical analysis of RTN was conducted locally on 100 devices of the small geometry W = 0.08 µm, L = 0.03 µm and reproduced on two different dies with two different sampling frequencies, which adds up to 400 samples per V g . An example of time-domain noise measurement using the addressable test structure and showing RTN is illustrated Fig. 11, with its associated amplitude ∆I d and trap mean capture and emission time constants τ c and τ e . Employing the processing of noise signals using smoothed histograms with kernel density estimation, as developed in [START_REF] Gauthier | Enhanced statistical detection of random telegraph noise in frequency and time domain[END_REF], an estimation of the distribution of active traps responsible for RTN was performed and is shown Fig. 12. The RTN-inducing trap distribution follows a Poisson law [START_REF] Nagumo | New analysis methods for comprehensive understanding of Random Telegraph Noise[END_REF], with the probability mass function given by:

f (k, λ) = λ k e -λ k! ( 3
)
where k is the number of RTN-inducing traps in our case, and λ is the expected rate of occurrence or mean. Poisson fits conducted on the distributions show a mean λ between 0.33 and 0.56 depending on the gate voltage, with a net increase for V g = 1 V, representing a proportion of RTN signatures ranging from 30% to 40% in the studied samples. This means that in the case of an analysis carried out globally with around 60 devices, while an estimation of the trap distribution is feasible, only a few RTN signatures can be found which is not sufficient to establish relevant distributions of trap parameters.

The addressable test structure thus allows the acquisition of a greater number of samples for the statistical analysis of RTN.

A cartography of the RTN-induced traps within the scribe for the two different dies at V g = 1 V is shown Fig. 13, where the local distribution of traps within the scribe seems random. Traps mean capture and emission time constants τ c and τ e were extracted and presented on a capture/emission time map Fig. 14, showing a single explicit cluster, meaning that the traps responsible for RTN are from the same type and nature [START_REF] Grasser | Analytic modeling of the bias temperature instability using capture/emission time maps[END_REF]. The distributions of the RTN amplitudes in terms of V t shifts, ∆V t = ∆I d /g m , are represented Fig. 15, showing possible log-normal distribution as mentioned in other studies [START_REF] Sonoda | Discrete Dopant Effects on Statistical Variation of Random Telegraph Signal Magnitude[END_REF]. The increase in ∆V t as V g increases could be explained by the general higher noise level in strong inversion associated with CMF, also discussed in section III-B. The intermediate V g shows a more scattered distribution, possibly due to the responsible traps being active in a wider gate voltage range, taking into account traps in common with the ones detected at V g = 0.3 V and V g = 1 V. The addressable test structure thereby allows a significant increase in the number of measured devices, consequently improving the statistical analysis of RTN, since its study needs numerous samples that can be scarce depending on the geometry or the technology.

IV. CONCLUSION

The viability of using an addressable array test structure for the statistical analysis of static mismatch, LFN and RTN has been demonstrated. Noise level and variability analysis performed on devices from the test structure were presented. Local and global variability analysis are made possible, showing a significant dominance of the impact of local variations on the variability. Local studies using the addressable test structure on a few dies can thus be conducted, allowing the increase in the number of samples and enabling a better understanding of the variations in static and noise parameters such as the number of traps and their characteristics. 571DPSOLWXGH∆V t P9 &XPXODWLYHSUREDELOLW\ V g 9 V g 9 V g 9

Fig. 1 .

 1 Fig. 1. Schematic representation of the addressable test structure.

Fig. 3 (Fig. 3 .

 33 Fig. 3(a) and Fig. 3(b) show the global I d -V g characteristics die-to-die spread of device #1 for the large and smaller

Fig. 4 .

 4 Fig. 4. Local (single die) transistor pairs mismatch standard deviation as a function of the inverse square root of device area.

Fig. 5 .

 5 Fig. 5. Transistor pairs mismatch standard deviation as a function of the inverse square root of device area for the single die local samples and the die-to-die global samples. The distribution of σ∆Vt for the smallest geometry is represented in the inset, with normal percentiles scale.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Input-referred gate voltage noise PSD S Vg of device #1 for all dies, of W = 0.5 µm, L = 0.5 µm transistors (a), and W = 0.08 µm, L = 0.03 µm transistors (b) for the same gate voltage Vg = 0.4 V.

Fig. 8 .*DWHYROWDJHV g 9 σFig. 9 .Fig. 10 .

 89910 Fig. 8. Median input-referred gate voltage noise PSD S Vg extracted at f = 10 Hz as a function of the drain current, for all measured devices. The global samples, local samples and the results from the independent standard test structure are compared to each other, showing no significant difference.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Example of time-domain drain current noise measurement showing a RTN signature, and its associated histogram.

Fig. 13 .

 13 Fig. 13. Cartography of the detected RTN traps within the columns of the test structure corresponding to the W = 0.08 µm, L = 0.03 µm devices, for two different dies, extracted at Vg = 1 V.

9 Fig. 14 .

 914 Fig. 14. Capture/emission time map extracted from the detected RTN signatures, for three different gate voltages.

Fig. 15 .

 15 Fig. 15. RTN amplitude ∆Vt distributions for three different gate voltages, with normal percentiles scale.

d DUT sense sense sense force force force

	Fig. 2. Circuitry surrounding each DUT of the test structure. Kelvin accesses
	are available for the drain, gate and source terminals. Each DUT terminal is
	selected through switches.
	addressing setup. Firstly, time and frequency domain noise
	measurements are done using the NOISYS equipment by
	Synergie Concept consisting of a low noise current amplifier,
	a precision voltmeter connected to the DC output of the
	amplifier, and a HP35670A spectrum analyzer connected to the
	AC output of the amplifier. The drain, gate, source, and bulk
	terminals of the DUT are connected to the bias units of the
	NOISYS equipment. Secondly, the addressing operation takes
	advantage of a switching matrix in which a voltage source