{'fraction_non_alphanumeric': 0.03046365005858394, 'fraction_numerical': 0.017463594264353066, 'mean_word_length': 4.754093097913323, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A theoretical approach for quantifying the impact of changes in effective population size and expression level on the rate of coding sequence evolution

Molecular sequences are shaped by selection, where the strength of selection relative to drift is determined by effective population size (Ne). Populations with high Ne are expected to undergo stronger purifying selection, and consequently to show a lower substitution rate for selected mutations relative to the substitution rate for neutral mutations (ω). However, computational models based on biophysics of protein stability have suggested that ω can also be independent of Ne, a result proven under general conditions. Together, the response of ω to changes in Ne depends on the specific mapping from sequence to fitness. Importantly, an increase in protein expression level has been found empirically to result in decrease of ω, an observation predicted by theoretical models assuming selection for protein stability. Here, we derive a theoretical approximation for the response of ω to changes in Ne and expression level, under an explicit genotype-phenotype-fitness map. The method is generally valid for additive traits and log-concave fitness functions. We applied these results to protein undergoing selection for their conformational stability and corroborate out findings with simulations under more complex models. We predict a weak response of ω to changes in either Ne or expression level, which are interchangeable. Based on empirical data, we propose that fitness based on the conformational stability may not be a sufficient mechanism to explain the empirically observed variation in ω across species. Other aspects of protein biophysics might be explored, such as protein-protein interactions, which can lead to a stronger response of ω to changes in Ne.

In this context, the aim of the present study is to characterize the dynamics and equilibrium response of 114 ω to changes in N e and expression level, and to relate this response to structural parameters of the model. 115 To this effect, we develop a general mathematical approach to derive a quantitative approximation of the 116 response of ω to changes in N e and expression level, in the context of a given genotype-phenotype-fitness map, 117 as depicted in figure 1. In the light of previously published empirical estimates from protein thermodynamics 118 and comparative genomics, we discuss the articulation between empirical data and our mechanistic model. 119 We also discuss some of the alternative biophysical mechanisms that could determine the selective landscape 120 on protein-coding sequences, and how they would modulate the response of ω to changes in N e and expression 121 level. The phenotype (x) is a real-valued function of the genotype (i.e. the amino-acid sequence), and is defined in our model as the fraction of destabilizing amino acids in the sequence. Fitness is a decreasing log-concave function of the phenotype, depending on structural parameters of the model. Once the relation from genotype to fitness is defined, the substitution process proceeds as presented in the middle panel. For a given effective population size N e , the evolutionary process results in an average value of the phenotype x * and an average substitution rate (relative to the neutral rate) ω. Averaging over time is equivalent to determining the statistical equilibrium, by ergodicity of the stochastic process. The slope of the scaling of the equilibrium ω as a function of log-N e defines the susceptibility χ, which is a function of the structural parameters defined by the phenotype-fitness map.
In the original biophysical model, protein stability is determined by the difference in free energy between 132 the folded and unfolded conformations, called ∆G and measured in kcal/mol. Technically, free energy is 133 computed based on the 3D conformation of the protein and using statistical potentials. As a result, the  Here, we approximate this model such that the (de-)stabilizing effect at a particular site, such as measured by the ∆∆G of the mutation, does not depend on other neighbouring residues, thus disregarding specific epistasis (Dasmeh et al., 2014). Instead, each site contributes independently and additively to ∆G. In addition, we assume that, for each site of the sequence, only one amino acid is stabilizing the protein. All 19 other amino acids are equally destabilizing. Each site bearing a destabilizing amino acid contributes an excess of ∆∆G > 0 (in kcal/mol) to the total ∆G. The smallest achievable value of ∆G, obtained when all amino acids of the sequence are stabilizing, is noted ∆G min < 0. In this model, the most succinct phenotype of a given genotype (i.e. sequence) is just the proportion of destabilizing amino acids in the sequence, defined as 0 ≤ x ≤ 1. Thus ∆G is a linear function of x: where n is the number of sites in the sequence.

138
For a given ∆G, thermodynamic equations allows one to derive the proportion of protein molecules that 139 are in the native (folded) conformation in the cytoplasm. This fraction is assumed to be a proxy for fitness, 140 motivated in part by the fact that a protein must be folded to perform its function. A slightly different model 141 will be considered below, in order to take into account protein expression level (see section 2.3).

142
Analytically, the fitness function is given by the Fermi Dirac distribution and is typically close to 1, leading to a first-order approximation (Goldstein, 2011): where W is the Wrightian fitness for a given phenotype and f is the Malthusian fitness (or log-fitness). Here,

143
∆G min and ∆∆G are defined as above, and the parameter β is 1.686 mol/kcal at 25°C (or 298.2K).

144
Of note, even though the phenotypic effect of a mutation at a given site does not depend on the amino-acids 145 that are present at other sites (i.e. the trait is additive), the fitness effect of a mutation still depends on other 146 sites (i.e. the log-fitness is not additive). As a result, the molecular evolutionary process is site-interdependent,  to discard the influence of the underlying mutation rate, is denoted ω. This relative rate can also be interpreted 152 as the mean fixation probability of mutations scaled by the fixation probability of neutral alleles p = 1/2N e , 153 the mean being weighted by the probability of occurrence of mutations in the population. As a result, an 154 ω < 1 indicates that mutations are negatively selected on average, and ω decreases with increasing strength 155 of purifying selection.

156
In this section we present an analytical approximate solution for the response of ω after a change in N e (in log space), as depicted in the right panel of figure 1. We call this response the susceptibility of ω to changes in N e , and denote it as χ: Deriving χ is done in two steps. First, we determine the mean phenotype at equilibrium, when evolutionary 157 forces of mutation, selection and genetic drift compensate each other. Subsequently, differential calculus is 158 used to compute the response of the equilibrium phenotype to a change in N e , which allows us to ultimately 159 derive an equation for χ. The main results of our derivation are given both in the general case of any and expression level on the rate of coding sequence evolution.
(log-concave) phenotype-fitness map, and in the specific case of the biophysical model introduced above. A 161 more detailed derivation is available in the supplementary materials.

162
For a given genotype, mutations can have various effects: they can increase or decrease the proportion of destabilizing amino acids, or do nothing if the mutation is between two destabilizing amino acids. To derive the probabilities of such events to occur, we also make the simplifying assumption that all transitions between amino acids are equiprobable. Altogether, any mutation in the sequence can then have a phenotypic effect of 0 or δx = 1/n, with probabilities of transitions equal to: In the extreme case of an optimal phenotype (x = 0), only destabilizing mutations are proposed. Moreover, 163 the probability to propose a stabilizing mutation (effect −δx), or a neutral mutation (effect 0), is proportional 164 to x. Conversely, the probability to propose a destabilizing mutation is equal to (1 − x). As a result, the 165 mutation bias is proportional to (1 − x)/x. This mutation bias fundamentally reflects a combinatorial effect, 166 due to the number of mutational opportunities available in either direction.

167
Second, we need to determine the strength of selection acting on mutations. Destabilizing mutations are selected against with a negative selection coefficient which can be approximated by: where f = ln(W ) is the log-fitness (or Malthusian fitness). Conversely, stabilizing mutations will be under 168 positive selection with opposite sign but same absolute value. It is important to realize that the selective 169 effect is dependent on x. Furthermore, because the fitness function is log-concave, the absolute value of s 170 increases with x.

171
Based on these expressions for the mutational and selective pressures, one can then study the trajectory followed by the evolutionary process. Starting from an optimal sequence, mostly destabilizing mutations will occur, some of which may reach fixation and accumulate until selection coefficients against new deleterious mutations is too strong, at which point the protein will reach a point of equilibrium called marginal stability (Taverna and Goldstein, 2002;Bloom et al., 2007). Most importantly, the probability of fixation of mutations is affected by genetic drift, and thus depends on the effective population size (N e ). At the equilibrium between mutation, selection and drift, the process fluctuates through the occurrence of advantageous and deleterious substitutions compensating each other. This equilibrium can be determined by expressing the constraint that the selection coefficient of substitutions is expected to be null on average (Goldstein, 2013).
Formally, and after simplification, the equilibrium phenotype denoted x * is given in the general case by: This equation cannot be solved explicitly for x * , but a qualitative intuition on the consequences of change 175 in N e to the equilibrium phenotype x * is given in figure 2. Intuitively, an increase in N e results in a more (right-hand side of eq. 10) increases exponentially with x where βn∆∆G is the exponential growth rate (yellow and green curves). When βn∆∆G is large, increasing N e by an order of magnitude (yellow dotted curves) very moderately impacts the equilibrium phenotype (small ∆x * ). In contrast, for small βn∆∆G (green curves), the equilibrium phenotype is more strongly impacted by a change in N e (large ∆x * ). Finally, response of x * to changes in N e reflects the response of ω since both are approximately equal (equation 11).
The results obtained thus far only relate the equilibrium phenotype (x * ) to N e . To capture how ω varies with N e , we also need to obtain an expression for ω as a function of x * . At equilibrium we can derive (supplementary materials) the expected substitution rate of mutations, and thus ω, which simply approximates to: (which are neutral) compose the largest proportion of proposed mutations having a substantial probability of 184 fixation (equation 6). In contrast, stabilizing mutations are rare, while destabilizing mutations have a low 185 probability of fixation. Since there is a fraction x * of sites already occupied by a destabilizing amino-acid, 186 these neutral substitutions occur at rate x * .

187
Combined together, these analytical approximations yield the susceptibility (equation 5) of ω to a change in N e : The two terms of the denominator correspond to the derivative of the mutational bias and the scaled selection coefficient, respectively. However, the mutational bias decreases weakly with x (blue curve on figure 2) while the strength of selection increases sharply with x (red and green curves). As a consequence, the derivative of the mutational bias is much lower than the derivative of the selection coefficient around the equilibrium point (i.e. the phenotype is nearly equimutable). The second term can therefore be ignored, which leads to a very compact equation for susceptibility χ in the general case: The susceptibility is thus equal to the inverse of the relative curvature, i.e. the ratio of the second to the first 188 derivatives, of the log-fitness function, taken at the equilibrium phenotype. Of note, this susceptibility is 189 strictly negative for decreasing log-concave fitness functions, asserting that ω is a decreasing function of N e .

190
In addition, the susceptibility itself is low in absolute value (i.e. ω responds more weakly) for strongly concave

194
In the specific case of the biophysical model, the susceptibility (χ) further simplifies to: meaning that ω is linearly decreasing with N e (in log scale) since χ is independent of x * , or, in other words,

207
Our general derivation can be directly applied to this case. For a given protein with expression level y and a cost A representing the selective cost per misfolded molecule (positive constant), the fitness and selection coefficient can be defined as follows: Under this model, the total selective cost of a destabilizing mutation is now directly proportional to the 208 total amount of misfolded proteins. This fitness function leads to the following expression for the mutation-209 selection-drift equilibrium: Importantly, in this equation, N e and y are confounded factors appearing only as a product. This means that increasing either N e or y leads to same change in equilibrium phenotype, and hence the same change in ω. In other words, the susceptibility of the response to changes in either N e or expression level is the same:   These simulation experiments demonstrate, first, that the relation between ω and log-N e is indeed linear, 237 at least in the range explored here, and that the slope of the linear regression matches the expected theoretical 238 value ( figure 3A). Secondly, we observe that the parameter ∆G min has virtually no effect on the slope of the 239 linear regression, as also expected theoretically ( figure 3B). Instead, decreasing ∆G min (to more negative 240 and expression level on the rate of coding sequence evolution. values) merely results in an overall increase in ω over the whole range of N e (i.e. has an impact on the intercept, 241 not on the slope of the relation). This is due to the fact that decreasing ∆G min shifts the equilibrium to 242 higher x * , since more destabilizing sites can then reach fixation before reaching the point of marginal stability.

243
Finally, we relaxed our assumption that each site of the sequence contributes independently to ∆G, by 244 taking into account the 3D structure of protein and using a statistical potential to estimate ∆G (supplementary closer to optimal and the purification of deleterious mutations is stronger. The reaction to a decrease in N e is A theoretical approach for quantifying the impact of changes in effective population size and expression level on the rate of coding sequence evolution. Overall, the susceptibility (χ) is a function of the structural parameters of the protein and takes a very 292 simple analytical form, being inversely proportional to the product of three terms: the sequence size, the 293 inverse temperature (β), and the average change in conformational energy of destabilizing mutations (∆∆G).

294
Quantitatively, this product can be several orders of magnitude greater than 1 in practice, such that the 295 susceptibility of ω, which is its inverse, is typically small. Previous studies using this model presented an 296 apparent lack of response of ω to changes in N e (Goldstein, 2013). We refine this result, by observing 297 that there is in fact a very subtle and weak relation, which requires extensive computation to be detected,  mutations is independent of ∆G is a necessary condition to observe independence between ω and N e (Cherry, 327 1998). In our model, the average ∆∆G of mutations at equilibrium depends on ∆G due to combinatorial 328 considerations, but this dependence is weaker than empirically observed, which also translates into a weaker 329 susceptibility of ω to changes in N e or expression level than empirically observed.

330
Thus, overall, the response of ω to either N e or expression level predicted by the biophysical model where W is the Wrightian fitness for a given phenotype and f is the Malthusian fitness (or log-fitness).

379
The waiting time before the next mutant invading the population, and the specific mutation involved where µ S,S is the mutation rate between S and S , determined by the underlying 4x4 nucleotide mutation rate Under the additive model for the free energy, the difference in free energy between folded and unfolded state is assumed to be given by: where 0 ≤ x (S) ≤ 1 is the distance of S to the optimal sequence (i.e. the fraction of sites occupied by a 388 destabilizing amino-acid). For each site of the sequence, the optimal amino acids are chosen randomly at 389 initialization, and the distance between the current amino acid and the optimal is scaled by the Grantham 390 amino-acid distance (Grantham, 1974). The Wrightian fitness is defined as the probability of our protein to where β is the inverse of the temperature (β = 1/kT ).

393
For simulations under a 3D model of protein conformations, we adapted the model developed in Goldstein 394 and Pollock (2017) to our C++ simulator (see supplementary materials).

395
For simulations under a site-independent fitness landscape, with site-specific fitness profiles, the protein 396 log-fitness is computed as the sum of amino-acid log-fitness coefficients along the sequence. In this model,