Exploiting the Adaptation Dynamics to Predict the Distribution of Beneficial Fitness Effects

Adaptation of asexual populations is driven by beneficial mutations and therefore the dynamics of this process, besides other factors, depends on the distribution of beneficial fitness effects. It is known that on uncorrelated fitness landscapes, this distribution can only be of three types: truncated, exponential and power law. We performed extensive stochastic simulations to study the adaptation dynamics on rugged fitness landscapes, and identified two quantities that can be used to distinguish the underlying distribution of beneficial fitness effects. The first quantity studied here is the fitness difference between successive mutations that spread in the population, which is found to decrease in the case of truncated distributions, remains nearly a constant for exponentially decaying distributions and increases when the fitness distribution decays as a power law. The second quantity of interest, namely, the rate of change of fitness with time also shows quantitatively different behaviour for different beneficial fitness distributions. The patterns displayed by the two aforementioned quantities are found to hold good for both low and high mutation rates. We discuss how these patterns can be exploited to determine the distribution of beneficial fitness effects in microbial experiments.


Introduction
Microbial populations have to constantly adapt in order to survive in a changing environment. For example, a bacterial population exposed to a new antibiotic must evolve in order to exist [1]. In asexual populations, this process of adaptation is driven only by rare beneficial mutations [2] which provide fitness advantage. Therefore, in order to survive in a new environment, enough beneficial mutations should be available and the beneficial mutations should confer sufficient fitness advantage. While the first factor depends on the mutation rate and population size, the second factor is determined by the underlying fitness distributions. Even though we have some understanding about the mutation rate of different microbial populations, the full fitness distribution is more complex and relatively little is known about it. However, for moderately adapted populations (i.e., fitness of the wild type is high enough), rare beneficial mutations which occur in the tail of the fitness distribution can be described by the extreme value theory (EVT) as proposed first by Gillespie [3]. The EVT states that the extreme tail of all distributions of uncorrelated random variables (fitness, in this case) can be of only three types. Depending on whether the tail of underlying fitness distribution is truncated or decaying faster than a power law or as a power law, the EVT distribution would belong to the Weibull or Gumbel or Fréchet domain, respectively [4]. All three EVT domains can be obtained from the generalized Pareto distribution given as where κ is the tuning parameter. One example from each of the three EVT domains is shown in Fig 1, which shows the distribution of beneficial effects p(f) with fitness f. The three types of EVT domains are classified according to the value of κ. Here negative κ belongs to the Weibull domain, while κ = 0 corresponds to the Gumbel domain and positive κ to the Fréchet domain. Interestingly, all three distribution of beneficial fitness effects(DBFEs) have been observed in experiments on microbial populations [5][6][7][8][9][10][11][12][13][14]. While the exponential distribution belonging to the Gumbel domain has been most commonly seen [5][6][7][8], in recent times, the distribution of beneficial mutations belonging to the Weibull [10,14] and Fréchet [11] domains have also been observed. Recent theoretical studies have shown analytically and numerically that qualitatively different patterns occur in the adaptation dynamics of populations in different EVT domains of DBFEs in a low mutation regime [15][16][17][18]. Specifically, it has been shown that fitness gain in a fixation event follows the pattern of diminishing returns in the Weibull domain, constant returns in the Gumbel domain and accelerating returns in the Fréchet domain, and thus indicates that this quantity can be used to predict the DBFE. These observations are restricted to strong selection-weak mutation (SSWM) regime in which the genetic variation in the population is minimal, that is, only one beneficial mutation is present in the population in the time interval between its appearance and fixation [7]. It is then natural to ask whether the relationship between adaptation dynamics and the DBFE mentioned above are robust for large populations, where there might be more than one beneficial mutation competing for dominance in the population. The main aim of our study is to address this question and to see if the fitness gain in a fixation event can be used for predicting the DBFE in a more general scenario.
Here, we are mainly concerned with the populations in which a large number of mutants are produced at every generation. Hence, more than one beneficial mutation is expected to be present at the same time [19][20][21][22][23]. In this case, the beneficial mutations will compete with each other as has been observed in different experimental populations [24][25][26][27]. In this high mutation regime, as a result of the competition among the beneficial mutations, the rate of adaptation slows down. Fitness advantage due to the mutations that get fixed is much higher, since the availability of more mutations results in allowing only the best (fittest) mutation to get fixed [28]. A clear comparison of the population fraction of new mutants appearing in a population for two mutation regimes is given in Fig 2. In Fig 2(a) we see that the population in the SSWM regime is more or less monomorphic with only one mutant present at a time in all the three EVT domains. However, in a high mutation regime, the population is polymorphic with more than one mutant produced in it at every generation as shown in Fig 2(b). In fact, a large amount of genetic variation is observed in the case of bounded distributions corresponding to κ < 0 in Eq (1) resulting in a strong competition between the beneficial mutants.
In this work, we have used Wright-Fisher dynamics to study the adaptation dynamics of an asexual population in high and low mutation regimes for the three EVT domains of DBFE. The main motivation of this study is to look for quantities which can be used to distinguish between DBFEs using the properties of adaptation dynamics as opposed to the direct measurements of DBFEs. Our most important and interesting result is concerned with the fitness difference between mutations that spread in a population. This quantity shows qualitatively different trends in three EVT domains and thus helps in distinguishing the DBFEs.
We have also studied another quantity which is the rate of change of fitness with time, and observed that this shows quantitatively different behaviour for different EVT domains of the DBFEs. Though some results for the rate of change of fitness are already known in the literature [29], we measured it for all the three cases (Weibull, Gumbel and Fréchet) and identified that this can be used to distinguish the DBFEs in both SSWM and high mutation regimes. In order to obtain a complete picture, a comparison of our study with the existing literature is given in Table 1 below.
We also measured quantities like the genetic variation and the number of mutations in the most populated sequence. All of these quantities are discussed in the Results section. We suggest that the distinct trends shown by the above mentioned quantities can be used to predict DBFEs from experimental studies on adaptation. The relevance of our work to experiments is also explored in the Discussion section.

Materials and Methods
We track the dynamics of a population of self-replicating (asexual), infinitely long binary sequences of fixed size using the standard Wright-Fisher process [21,28]. In our work, the population size is held constant at N = 10 4 , unless specified otherwise and the total mutation probability (beneficial and deleterious) per sequence is given by μ. Every occupied sequence is counted as a class and is labeled when it arises in the population. Initially, the whole population is in class 1 whose fitness is fixed and specified in every simulation run. We have used the term leader to refer to the class whose normalised probability of reproduction (product of population fraction and fitness) is greater than half. In that case, clearly class 1 is the initial leader since the whole population is localized there. At every time step, out of N sequences, m t are chosen from a binomial distribution with mean Nμ as mutants. Every mutant produced increases the number of classes in the population by one, and with time, the mutants may produce their own set of further mutants. The population fraction of each class may grow or go extinct, as can be observed in Fig 2. At any time t, the number of classes present in the population is given by N c ðtÞ, and the population size and fitness of each class, i, where 1 i N c , is  Note that though the fitness of the class is the same as long as it persists in the population, its size may vary at every time step, thus changing its probability of reproduction as given by Eq (2). Different classes are populated in the next time step based on the multinomial distribution where t 0 = t + 1. The above equation is subject to the constraint S N c ðtÞ j¼1 nðj; t 0 Þ ¼ N. In our simulations, we implement Eq (3) along with the above constraint by converting Eq (3) We set the population size of the last class as nðN c ðtÞ; andÑ ðiÞ ¼ N À S iÀ1 j¼1 nðj; tÞ. At every time step, once the classes are populated based on the algorithm described above, m t sequences are chosen as mutants based on the binomial distribution with mean Nμ. Every new mutant class that appears in the population reduces the population size of the class in which it arose by one. In our work, we have varied μ to access both the SSWM (low mutation) and the high mutation regime. In our simulations unless specified otherwise, Nμ = 0.01 in low (SSWM) and Nμ = 50 in high mutation regimes.
A new class is assigned to each mutant and its fitness is chosen from a generalized Pareto distribution [4] given in Eq (1). The advantage of using Eq (1) is that we can access all three EVT domains of DBFE by changing κ. The distributions whose κ < 0 belong to the Weibull domain, while κ = 0 belong to the Gumbel domain, and κ > 0 belong to the Fréchet domain, respectively. The frequency distribution of beneficial effects p(f) for various values of κ is shown in Fig 1. The upper bound u for the distributions chosen from Eq (1) is infinity when κ ! 0 and equals −1/κ for κ < 0. In this work, the fitness of the mutants is independently chosen from Eq (1) thus making the fitness of the mutant, F m an uncorrelated variable, which may be greater or smaller than the parent fitness, F p . We have analyzed the results to see how they vary between the three EVT domains and different mutation rates.
In the allocation of the fitness to any mutant, our work differs from the other works on clonal interference [21,28] wherein the fitness of the mutant is hiked above the parent fitness by the selection coefficients (s) which may be held constant or chosen from a distribution as F m = (1 + s)F p . Unlike the model we have used in this work (as explained above), in this case, there is a strong correlation between the mutant fitness F m and the parent fitness F p . In those cases, the mutant fitness is always greater than the parent fitness and on an average, a double or higher mutant is fitter than a single mutant. This is in contrast with our work since in ours, as the fitness of the parent increases, the number of better mutants available decreases thus producing different patterns for the fitness increment in each EVT domain.
In our model, whenever a mutant class goes extinct, the classes below it are moved up and the number of classes in the population is reduced by one. The normalised probability of reproduction given in Eq (2) of a class exceeding half corresponds to a leader change. The new leader determined now belongs to the class whose normalised probability exceeded half. We have also explored other criteria for defining the leader as the most populated class and find that our main results are robust with respect to the change in criteria (data not shown).
Every change of a leader is counted as a step. In the high mutation regime, the population is spread over many sequences and a sequence can produce two or more mutants, each of which may become leaders at different time steps. However, in the SSWM regime, the whole population is localised at a single sequence with a fixed fitness and can only move to a different sequence with higher fitness one mutation away. Thus every new leader arises from the previous leader, as can be observed in Fig 2(a). When a better sequence appearing in the population does not get lost due to genetic drift, it quickly gets fixed. Further mutations that may lead to future leaders appear in this genetic background. The change in the fitness of the population is the same as the change in the fitness of the leader. In this case, every move of the population (leader) from one sequence to another is termed as a step in the adaptive walk [30][31][32][33], whereas in the high mutation regime, the population is polymorphic and as seen from Fig 2(b) the leader change is not obvious.
Various quantities like the difference in fitness between successive leaders and the average number of mutations in the leader are averaged only over the walks that take the step. Other quantities like the number of classes present at any point in time and the rate of change of fitness are averaged over all time steps in that simulation run.
In this paper, the total number of iterations is 10 5 in every simulation run and the dynamics are tracked for a finite time limit of 10 4 generations, which we shall refer to as t max . In this time span, the maximum fitness value, f max that arises in the population can be calculated as where u is the upper limit of the fitness distribution equalling (-1/κ) for bounded distributions and infinity for the unbounded ones [4]. From the above integral, we get Results

The number of classes in the population
For a population that is fixed in size, the number of classes in the population is expected to increase with the mutation rate. The average genetic variation, which is defined here as the average number of classes (N c ) present in the population is shown in Fig 3 for all three DBFE domains. The top and bottom panels of the figure show the data corresponding to the high and low mutation regimes respectively. In both mutation regimes, we see that the average number of classes increase during the initial time steps and decrease at later times when the classes with lower fitness are eliminated by the fitter ones. The maximum number of classes existing in the population for the first case, as shown in Fig 3(a), does not belong to the lowest initial fitness, but to a slightly higher initial fitness. This could be because when the initial fitness is low, its class is quickly replaced by a fitter mutant and all further mutants that arise on this new background must compete with this fitter class.
In the low mutation regime, the population is localized at a single sequence for most of the time and produces Nμ mutants at every time step. Hence, in this case, the average number of classes approach a constant Nμ + 1 at large times as can be seen in the bottom panels of Fig 3. These panels also indicate that the value of this constant increases with decreasing κ. This is because in the case of bounded distributions with κ < 0, the fitness of a beneficial mutant produced is expected to be closer to the parent fitness. In other words, mutations are nearly neutral and thus it takes a longer time to take over the population as shown in Fig 2(a). This results in a larger number of mutants in the Weibull domain, which can be observed in the bottom panel of Fig 3(a). We can clearly see from the top panels of Fig 3 that number of classes increases with decreasing κ even in a high mutation regime. Moreover, the average number of classes present at a time is much higher in this regime. This makes sense because the fitness of the classes belonging to κ = −1 cannot be very different from each other (can take on values between 0 and 1), which makes it possible for many of them to exist in the population. The maximum fitness of the classes belonging to κ = 1/4 distribution will on an average be much higher than all others (since the distribution is unbounded with a fat tail), thus out-competing the others in the population.

Number of mutations in the leader
In the low mutation regime, the average number of mutations in the leader is expected to be very close to the step number since the genetic variation in the population is low and any mutation that escapes drift quickly takes over the population [3]. We verify this point via simulations as depicted in Fig 4. We find that the mutation number equals the step number in all the three EVT domains of the DBFE in the low mutation regime during the initial steps. However in the high mutation regime, the number of mutations in the leader of any step differs between the three DBFE domains. When the mutation rate is increased, the genetic variation of the population and the significance of clonal interference also increases. In the high mutation regime, the number of mutations in the leader is found to be less than the step number in all three DBFE domains. This is because there is a chance that different mutants originating from the same parent class can become the leader of the population at different times. This decrease from the step number is the minimum for the fat-tailed distributions and maximum for the truncated ones, as shown in Fig 4. This result is consistent with the number of classes present in the population as discussed in the previous section. In the Fréchet domain, since the clonal interference is minimal, it is most likely that a mutant originating from the present leader will become the next one. In the Weibull domain, due to the large number of classes present in the population, mutants originating from the same class can become leaders at different time points.

Fitness and fitness difference
From our simulations, we find that the average fitness of the first mutant fixed in the population, f 1 increases linearly with initial fitness, f 0 for all κ in the low mutation regime and for κ 6 ¼ 0 in the high mutation regime. So we can write where the coefficients a ðNmÞ k and b ðNmÞ k are constants. In the low mutation regime, where the population for most times is monomorphic, the adaptive walk model has been used to analytically obtain the fitness at the first step, f 1 as [15,16] where the transition probability In this model, from Eq (9), the coefficient a ðNm(1Þ k was obtained as 0.33, 1.0 and 1.6 for κ = −1, 0, and 1/4, respectively. The corresponding b ðNm(1Þ k for the aforementioned κ were 0.66, 2.0 and 1.89 [16]. In the high mutation regime where the adaptive walk model is not applicable, we obtained the values for the coefficients in Eq (8) numerically. We find that for large f 0 , a ð50Þ k equals 0.004 and 1.5 and b ð50Þ k equals 0.99 and 9.1 for κ = −1 and 1/4 respectively. The interesting result from our work is that, irrespective of the number of mutants produced in the population, the difference Df step ¼ f 1 À f 0 between the fitness of the first step and the initial fitness displays different qualitative trends: it increases for positive κ, approaches a constant when κ = 0 and decreases for negative κ as shown in Fig 5 and S1 Fig. We can better understand these increasing and decreasing trends by the following heuristic argument. In both the low and high mutation regimes, for large f 0 , the fitness at the first step f 1 increases linearly with the initial fitness is given in Eq (8). Therefore, we can write the selection coefficient defined as the relative fitness difference at the first step as In an adapting population, since the fitness of the first step is greater than the initial fitness, the selection coefficient is always positive. As the fitness distributions belonging to the Fréchet domain are unbounded with fat tails, high f 0 values can be considered. In this case, the second term on the right hand side (RHS) of Eq (11) can be ignored and we can write s % ða ðNmÞ k À 1Þ > 0. Thus for κ > 0, since a ðNmÞ k > 1 it follows that the fitness difference at the first step increases with f 0 . On the other hand, since the distribution belonging to the Weibull domain is truncated, we can invoke the following inequality to explain the decrease in fitness difference with increasing f 0 : where u is the upper limit of the fitness distribution. With increasing f 0 , the RHS of the above equation decreases showing that as the initial fitness increases, f 1 À f 0 has to necessarily decrease. Thus, the qualitative trends discussed above appear to be determined by the behaviour of the tail (bounded/unbounded), and not by the details of the model. Further, it is interesting to note that while the data points for the exponentially decaying distribution (κ = 0) increase and seem to be approaching a constant in the low mutation regime, the data in the high mutation regime seems to be reducing to approach the same constant. Our  (9), while the dotted lines are the approximate results that can be obtained for the results when the initial fitness is high in the low mutation regime. The broken lines for κ 6 ¼ 0 are lines of best fit as mentioned in the text. The broken line for κ ! 0 is used for connecting the data points. The inset shows the fitness difference at the first step as a comparative measure of the fitness difference obtained at the first step when f 0 = 0. Here, the lines are used for connecting the data points. simulation results shown in Fig 5 not only match the predicted theoretical values and validate the claim of different qualitative trends in each EVT domain in the SSWM regime, but also show that the trends hold irrespective of the number of mutants produced in the population. This result suggests that the qualitatively different trends of the fitness difference (increasing, constant and decreasing with initial fitness in the Fréchet, Gumbel and Weibull domains, respectively), can be used to distinguish between the EVT domains in a more general scenario.
Though the fitness difference at the first step is greater in the high mutation regime, when compared with the results in the low mutation regime, when we look at the fitness difference at the first step scaled by the fitness difference obtained when the initial fitness is zero (insets of Fig 5), we see that this increase is slower in the high mutation regime compared to the results obtained in the low mutation regime. This indicates that as the mutation rate increases, though the number of mutants accessed is higher, the difference in fitness compared to a lower initial fitness is not proportionally higher and is in fact lower for all the fitness distributions.

Rate of change of fitness with time
Besides the fitness increment at a fixed event of leader change, we also measured the fitness as a function of time as shown in Fig 6. We observed that even though the fitness increases with time in all the three EVT domains, the rate at which the fitness increases depends strongly on the DBFE. This rate has an initial fast transient phase, after which it slows down.
The initial transient phase is strongly dependent on the initial condition as well as the mutation rate as shown in S2 Fig. The increase in fitness is fastest for the lowest initial condition, but it approaches the same fitness value as in the case of higher initial fitness in few generations. The time taken for populations of different initial fitness to reach the same fitness value depends on the mutation rate: for Nμ ) 1, it takes about 20 generations, whereas for Nμ ( 1,  Adaptation Dynamics Predicts the DBFE it is approximately 200 generations. Even after this transient phase, the rate of increase in average fitness ( F ðtÞ) with time depends on the mutation rate as shown in Fig 6. This is because of the fact that when a large number of mutations is available at the same time, a highly fit mutant can invade the population and give a large fitness increment. Therefore, the fitness of a highly fit mutant sequence would be greater in the high mutation regime compared to the one in a low mutation regime. The maximum fitness value reached in 9000 generations, in the case of Fréchet distribution, is about 10 times more for the high mutation regime, which is consistent with the expectation from Eq (7). Even beyond this point we noticed that the fitness is still increasing. In the same way, the Gumbel distribution also shows a significant increase in maximum fitness reached in the high mutation regime as compared to the SSWM regime (about 4 times). Here also we found that the fitness is still increasing beyond the time point till which we tracked the dynamics. The bounded distribution (Weibull) reaches near the upper bound in SSWM and evolves slowly. However, fitness reaches a fitness plateau in the high mutation regime and rate of adaptation becomes zero as can be seen in Fig 6. From this, we observe that the rate of change of fitness strongly depends on the properties of the underlying DBFE, which suggests that looking at this quantity can help us in distinguishing the DBFEs. Hence, we measured the fitness increment defined as at each step. The D F ðtÞ initially increases, then slowly decreases and settles down to a zero as shown in Fig 7. If we denote this function as where A is a constant and the exponent α can be used to distinguish the DBFE, since, as explained below, exponent α is found to be greater (smaller) than one in the Weibull (Fréchet) domain, but is close to one in the Gumbel domain.
In the SSWM regime, from Fig 7(a), we can see that each type of DBFE considered shows a different rate of decay. The Weibull domain has a faster decay with α = 1.86, the Gumbel domain has α % 1 [29] and the Fréchet domain α = 0.66 [29]. We observed that the same trend is robust in a high mutation rate regime as well, where α values are slightly larger in all cases. In this regime also α = 2.02, 1 and 0.76 for the Weibull, Gumbel and Fréchet domains, respectively as shown in Fig 7(b). In the high mutation regime, in the case of Weibull distributions, fitness reaches a plateau in few generations, after which its rate of change goes to zero as observed in Fig 7(b). The theoretical prediction for fitness at every time step for the unbounded distributions belonging to the Gumbel and Frèchet domains was obtained by Park and Krug [29] in the low mutation regime. The comparison of our simulation data with these predictions shows a very good agreement in the Gumbel domain and in the Fréchet domain (up to a constant). In this work, we have also considered the bounded distribution and observed that its rate of decrease is faster with an exponent greater than one, which was not considered in the previous studies. We observed that even in a high mutation regime, the exponent α shows the same behaviour. In this regime, the rate of change of fitness has been calculated only for exponential distribution belonging to the Gumbel domain [29] and their prediction matches with our data. In this work, we have obtained a complete picture by studying the rate of change of fitness numerically for the other two EVT domains as well.
Thus, the second main finding from our study is that in all DBFEs, the fitness difference at each time step decreases with time as given by Eq (14) and we can distinguish between the three EVT domains of DBFEs by looking at the exponent α.

Discussion
The main purpose of our work is to determine the quantities that can be used to distinguish the different extreme value domains of DBFE. Previous studies [16,18] have found that in an adapting population, the fitness gain at each fixation event shows qualitatively different trends in the three DBFE domains when the number of mutants produced in the population is much less than one at every generation (Nμ ( 1). The focus of this work is to explore the parameter regime in which the number of mutants produced is much above one (Nμ ) 1). When the mutation rate is high, the population becomes polymorphic and the better mutants existing in the population compete with each other. From our study, we have observed that the qualitative trends found for fitness difference when a new mutation establishes in the low mutation regime hold irrespective of the number of mutants produced. Thus, this study suggests that the  In each case the data is fitted with the theoretically expected function given in Eq (14), except for the exponential distribution for which we used the theoretical prediction by Park and Krug [29]. In all cases, the population starts with the same initial fitness f 0 = 0.5. fitness difference between successive mutations that spread in the population is a very important and robust quantity that can be used to predict the DBFEs in a more general scenario.
From our simulations, we see that as the initial fitness is increased, the fitness difference at the first step given by Df step reduces, approaches a constant, or increases with the initial fitness in the Weibull, Gumbel and Fréchet domains, respectively. We can understand these trends by a heuristic reasoning as discussed in detail in the Results section. This argument explains the increase in Df step with f 0 for an unbounded power law distribution and shows that the trends are determined by the behaviour of the tail (bounded/unbounded), and not by the details of the model.
Another important measure in understanding the dynamics of adaptation is the rate at which it occurs. Most of the previous studies which measured the adaptation rate have only considered exponentially distributed fitness distributions [20-22, 28, 34]. A previous study by Park and Krug [29] also considered DBFEs belonging to the Fréchet domain, but only in the SSWM regime (see Table 1). In this work, we have extended the previous studies by numerically measuring the rate of change of fitness for bounded distributions as well. We have measured the rate of change of fitness in all the three EVT domains of the DBFE in both low and high mutation regimes. We observed that in all the cases, the rate of change of fitness decreases with time as *t −α , where α > 1 for Weibull, α % 1 for Gumbel [29] and α < 1 for Fréchet domains [29].
Experimentally, the distribution of beneficial fitness effects can be inferred by two methods. In the first method, mutations are introduced in the wild type sequence and those that confer a fitness advantage are separated and their distribution of fitness effects are determined. In this method, DBFE belonging to all the EVT domains have been observed [5][6][7][8][9][10][11][12][13][14]. In contrast, here we focus on learning about DBFE via adaptation dynamics. Though many works have tracked the dynamics of the population during adaptation [7,[35][36][37][38], in most of them only the selection coefficient of the mutant fixed was measured. In our study, we have observed that the selection coefficient as given by Eq (11) always decreases, with the increasing initial fitness or increasing steps as shown in S3 Fig. Hence, this quantity is not useful to distinguish between the EVT domains. However, from our study we observe that the fitness difference between steps shows different patterns depending on the EVT domain of the DBFEs in both the high and low mutation regimes and can be used to distinguish between the EVT domains.
In this work, we have numerically shown that the fitness returns in each EVT domain is very robust and holds good even when the number of mutations produced is large (Nμ ) 1). Fitness difference can be measured in experiments, for example as in [5]. We suggest that experiments can predict the EVT domain of DBFE by measuring the fitness difference between successive mutations fixed in the population or even from the fitness of the first mutation, when the initial fitness is varied. However, currently experimental studies that measure both fitness and DBFE in the same study are not available, but it is highly desirable to have such studies to test our predictions.