Théo Roncalli
email: theo.roncalli@universite-paris-saclay.fr

Loïc Paulevé
email: loic.pauleve@labri.fr

Variable-Depth Simulation of Most Permissive Boolean Networks

In systems biology, Boolean networks (BNs) aim at modeling the qualitative dynamics of quantitative biological systems. Contrary to their (a)synchronous interpretations, the Most Permissive (MP) interpretation guarantees capturing all the trajectories of any quantitative system compatible with the BN, without additional parameters. Notably, the MP mode has the ability to capture transitions related to the heterogeneity of time scales and concentration scales in the abstracted quantitative system and which are not captured by asynchronous modes. So far, the analysis of MPBNs has focused on Boolean dynamical properties, such as the existence of particular trajectories or attractors. This paper addresses the sampling of trajectories from MPBNs in order to quantify the propensities of attractors reachable from a given initial BN configuration. The computation of MP transitions from a configuration is performed by iteratively discovering possible state changes. The number of iterations is referred to as the permissive depth, where the first depth corresponds to the asynchronous transitions. This permissive depth reflects the potential concentration and time scales heterogeneity along the abstracted quantitative process. The simulation of MPBNs is illustrated on several models from the literature, on which the depth parametrization can help to assess the robustness of predictions on attractor propensities changes triggered by model perturbations.

Introduction

Boolean networks (BNs) have been employed to model the temporal evolution of gene expression and protein activities in biological systems [START_REF] Kauffman | Metabolic stability and epigenesis in randomly connected nets[END_REF][START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF][START_REF] Zañudo | Cell line-specific network models of ER+ breast cancer identify potential PI3kα inhibitor resistance mechanisms and drug combinations[END_REF][START_REF] Montagud | Patient-specific boolean models of signalling networks guide personalised treatments[END_REF]. A BN is composed of a finite set of components having two states, either 0 (false) or 1 (true). The BN then specifies in which contexts the components can change state, such as "component 3 can switch to state 1 if and only if component 1 is in state 0 and component 2 is in state 1; in all other contexts, component 3 can only switch to 0". These rules can be given as one Boolean function per component, associating each possible configuration (which associates each component to a state) to a Boolean value: in our example, the function of component 3 is f 3 : {0, 1} n → {0, 1} with f 3 (x) = ¬x 1 ∧ x 2 , where ¬ and ∧ denote the logical negation and conjunction, respectively.

The execution of a BN relies on a given update mode which specifies how are computed the evolution of the states of components. In the systems biology literature, two main update modes are widely employed: the synchronous (or parallel ) and fully-asynchronous. In synchronous, each component is updated simultaneously: the configuration x = (x 1 , x 2 , x 3 ) evolves in one step to (f 1 (x), f 2 (x), f 3 (x)). In fully-asynchronous, only one component can be updated at a time, possibly leading to non-deterministic transitions: the configuration x can evolve either to (f 1 (x), x 2 , x 3 ), (x 1 , f 2 (x), x 3 ), or (x 1 , x 2 , f 3 (x)). Whatever the chosen update mode, as there is a finite number of configurations, an execution eventually reaches a set of limit configurations which will be visited infinitely often. These sets of configurations are called attractors: each execution eventually reaches and stays within one attractor. Attractors are prominent dynamical features when modeling biological processes: they represent stable behaviors and are usually associated to cellular phenotypes.

Computing the attractors reachable from a given initial configuration is at the core of many studies of biological processes with BNs [START_REF] Calzone | Mathematical modelling of cell-fate decision in response to death receptor engagement[END_REF][START_REF] Cohen | Mathematical modelling of molecular pathways enabling tumour cell invasion and migration[END_REF][START_REF] Remy | A modeling approach to explain mutually exclusive and co-occurring genetic alterations in bladder tumorigenesis[END_REF][START_REF] Mendes | Estimating attractor reachability in asynchronous logical models[END_REF][START_REF] Montagud | Patient-specific boolean models of signalling networks guide personalised treatments[END_REF]. These studies typically involve comparing the effect of a network mutation (forcing some components to have fixed value) on the sets of reachable attractors and their propensities. This later notion is often related to the number of paths leading from the initial configuration to each attractor: under some mutations, the same set of attractors may still be reachable, but the proportion of trajectories leading to them may substantially differ. This motivated the development of simulation algorithms for BNs in order to sample trajectories and quantify the propensities to reach attractors. These methods replace the non-determinism of asynchronous transitions by a probabilistic choice [START_REF] Mendes | Estimating attractor reachability in asynchronous logical models[END_REF][START_REF] Stoll | Continuous time boolean modeling for biological signaling: Application of gillespie algorithm[END_REF].

Whenever employed as qualitative models of quantitative systems, dynamics of BNs aim at giving a coarse-grained view of system dynamics without requiring numerous quantitative parameters. However, the Boolean (a)synchronous modes are not correct abstractions of quantitative dynamics [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF]: they lead at the same time to predict spurious transitions and, importantly, preclude transitions which are actually possible when considering delays for instance. Let us illustrate this with the following BN, denoted (A) in the rest of the text:

f 1 (x) = 1 f 2 (x) = x 1 f 3 (x) = (¬x 1 ∧ x 2 ) ∨ x 3 (A)
From the initial configuration where all components are 0, we write 000, there is only one possible transition: the activation of 1, leading to the configuration 100. Then, again, only one transition is possible, activating 2, thus leading to 110. There, the execution will stay infinitely on this configuration: no other state changes are possible: from 000 it is impossible to eventually activate 3, and {110} is the only reachable attractor. However, it is known that this system can actually activate 3 for a range of kinetics, as observed experimentally [START_REF] Schaerli | A unified design space of synthetic stripe-forming networks[END_REF], and easily captured with quantitative models [START_REF] Ishihara | Cross talking of network motifs in gene regulation that generates temporal pulses and spatial stripes[END_REF][START_REF] Rodrigo | Structural discrimination of robustness in transcriptional feedforward loops for pattern formation[END_REF]. Indeed, consider that 1 has actually several activation levels: one intermediate 1 /2 which is sufficient to activate 2 but not enough sufficient to inhibit 3: when in 1 /200, 2 can be activated, going to 1 /210, and then 3 can change state, going to 1 /211. When 1 becomes fully active (111), 3 will self-maintain its activation. Thus, a correct Boolean analysis of the BN f above should conclude that from 000, two attractors are reachable: {110}, when 1 goes rapidly to its maximum level, and {111} when 3 had time to activate before the full activation of 1. This example shows the limit of (a)synchronous interpretations of BNs when used as abstraction of quantitative systems: they enforce that the Boolean 0 matches with the quantitative 0, and the Boolean 1 matches with the quantitative non-zero (> 0), making impossible to capture transitions happening at different activation levels or different time scales.

The Most Permissive (MP) update mode of BNs [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF][START_REF] Paulevé | Non-deterministic Updates of Boolean Networks[END_REF] is a recently-introduced execution paradigm which enables capturing dynamics precluded by the asynchronous ones. The main idea behind the MP update mode is to systematically consider a potential delay when a component changes state, and consider any additional transitions that could occur if the changing component is in an intermediate state. It can be modeled as additional dynamic states "increase" ( ) and "decrease" ( ): when a component can be activated, it will first go through the "increase" state where it can be interpreted as either 0 or 1 by the other components, until eventually reaching the Boolean 1 state. With the previous example, starting from 000, the first component is put as increasing, thus going to the MP configuration 00. In this configuration, 2 can be activated because can be interpreted as 1, leading to 0. Then, 3 can be activated, as it can interpret the dynamic state of 1 as the Boolean 0, and 2 as the Boolean 1. This model the fact that the component 1 is not high enough for inhibiting 3, while 2 is high enough to activate it. Thus, in this example, there is a trajectory from 000 to 111, i.e., a configuration where all the components are active. As 111 is a fixed point of f , the MP analysis would thus conclude that two attractors are reachable from 000: {110} and {111}.

The MP update mode brings a formal abstraction property to BN dynamics with respect to quantitative models, without requiring additional parameters: essentially, MPBNs capture any behavior that is achievable by any quantitative model being compatible with the logic of the BN. This include, for instance, models which result from introducing quantitative parameters, such as transition speed and interaction thresholds. We give here a brief informal overview of the property (see [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF] for details): let us consider multivalued networks (MNs) of dimension n where components can have values in M = {0, . . . , m} for some m ∈ N >0 . A MN can be considered as map from multivalued configurations to the derivative of the value of the components, i.e., of the form F : M n → {-1, 0, 1}. A multivalued configuration z ∈ M n can be binarized by associating components with value 0 to the Boolean state 0, components with value m to the Boolean state 1, and other components to any Boolean state. A MN F is a refinement of a BN f whenever for any multivalued configuration z and for each component i, if component i increases (resp. decreases), i.e., F i (z) > 0 (resp. F i (z) < 0), there exists a binarization of this configuration such that f i is evaluated to 1 (resp. to 0). Then, for any pair of Boolean configurations x and y, if there exists a trajectory from m • x to m • y in the asynchronous dynamics of the MN, there necessarily exists a trajectory from x to y in the MP dynamics of the BN. Moreover, to any transition computed according the MP update mode, there is refinement of the BN which realizes it. One of the major consequences of the abstraction property is that if a configuration is not reachable from another with the MP update mode, then no quantitative models being compatible with the BN can produce the trajectory. In addition, the MP mode has a lower computational complexity for computing reachability and attractor properties, enabling formal analysis of genome-scale BNs [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF][START_REF] Paulevé | VLBNs -Very Large Boolean Networks[END_REF].

However, the simulation of MPBNs, i.e., the sampling of transitions following the MP update mode, has not been addressed so far. Thus, besides computing Boolean properties, such as the existence or absence of reachable attractors, there is no algorithm nor tools to approach the effect of a mutation on the propensities to reach attractors, as we mentioned above with asynchronous BNs.

In this paper, we present a first algorithm for sampling trajectories of BNs with the MP update mode, subject to additional simulation parameters for assigning probabilities to the transitions. The MP transitions enabled from a single configuration are computed iteratively, and we refer to the number of times this iteration is performed as the depth of the MP computation. At depth 1, the transitions match with the asynchronous update mode, but further depths bring additional behaviors. These iterations capture possibly different time scales: while component 1 is changing (depth 1), 2 can change (depth 2); then while 2 is changing, 3 can change (depth 3), etc. In our simulation algorithm, the probability of a transition can be affected by its depth and the number of components it changes simultaneously. Thus, similarly to [START_REF] Mendes | Estimating attractor reachability in asynchronous logical models[END_REF][START_REF] Stoll | Continuous time boolean modeling for biological signaling: Application of gillespie algorithm[END_REF] with the fully-asynchronous mode, our sampling of trajectories can be assimilated to a random walk in the MP dynamics, where the probability of transitions can be tuned with the simulation parameters. As we will show on case studies, the MP interpretation can lead to drastic changes in the predicted probabilities of reachable attractors, enabling assessing the robustness of prediction to the heterogeneity of time scales and concentration scales in the quantitative system captured by the discrete MP dynamics. Nevertheless, the simulation parameters and derived transition probabilities are empirical and cannot be formally related to a putative abstracted system.

Background

Boolean networks and dynamics

A Boolean network (BN) of dimension n is a function f : B n → B n , with B = {0, 1}. For each i ∈ {1, . . . , n}, f i : B n → B is the local function of compo- nent i.
The Boolean vectors x ∈ B n are the configurations of f , where for each i ∈ {1, . . . , n}, x i is the state of i. Given two configurations x, y ∈ B n , the components having a different state are noted ∆(x, y) = {i ∈ {1, . . . , n} | x i = y i }.

The influence graph G(f ) of a BN f is a signed digraph whose nodes are the components and edges mark the dependencies between them in the local functions: for all i, j ∈ {1, . . . , n}, there is an edge i s -→ j in G(f ) with s ∈ {-1, +1} if and only if there exists a configuration x ∈ B n with x i = 0 such that s = f j (x 1 , . . . , x i-1 , 1, x i+1 , . . . , x n ) -f j (x): the sole increasing of component i causes f j to increase (s = +1) or decrease (s = -1). Note there may exist i and j so that both i +1 --→ j and i -1

--→ j are in G(f ), for instance with f j (x) = x i xor x k . A BN f is locally monotone whenever there is no i, j ∈ {1, . . . , n} such that both i +1 --→ j and i -1

--→ j are in G(f ), i.e., each of its local function is unate. Computing G(f ) is a DP-complete problem (both in NP and coNP) [START_REF] Crama | Boolean Functions[END_REF]. For the example BN

f of Eq.(A), G(f ) = {1 +1 --→ 2, 1 -1 --→ 3, 2 +1 --→ 3, 3 +1 --→ 3}: it is locally monotone.
An update mode µ of f specifies a binary transition relation between configurations → µ ⊆ B n × B n . Classical update modes include the synchronous (or parallel ) update mode where x → s y iff x = y and y = f (x); the fully-asynchronous update mode where x → a y iff x and y differ on only one component i and y i = f i (x); and the (general) asynchronous update mode where x → g y iff x = y and for each component i ∈ ∆(x, y), y i = f i (x).

Given an update mode µ, a configuration y ∈ B n is reachable from a configuration x ∈ B n , noted x → * µ y, if and only if either x = y or there exists a sequence of transitions x → µ • • • → µ y. A non-empty set of configurations A ⊆ B n is an attractor if and only if for each pair of configurations x, y ∈ A, y is reachable from x, and there is no configuration z ∈ B n \ A that is reachable by a configuration in A. Remark that attractors are the bottom strongly connected components of the digraph (B n , → µ ). Whenever A is a singleton configuration {x}, it is said to be a fixed point of the dynamics. Otherwise, A is a cyclic attractor. In the case of (a)synchronous and MP update modes, the fixed points of the dynamics match exactly with the fixed points of f , i.e., the configurations x ∈ B n such that f (x) = x. Finally, the strong basin of an attractor A is the set of configurations that can reach A and no other distinct attractor.

Fig. 1 shows the transitions computed with the above defined update modes with the BN f defined as follows:

f 1 (x) = x 1 ∧ ¬x 3 f 2 (x) = x 1 f 3 (x) = ¬x 1 (B)
This model has two attractors, being fixed points 001 and 110. With the syn-chronous mode, only 001 is reachable from 111, whereas both are reachable with the asynchronous modes.

Sub-hypercubes and closures

The MP update mode and the algorithm presented in this paper gravitate around the notion of sub-hypercube of B n and their partial closure by f . A Boolean sub-hypercube of dimension n is specified by a vector in {0, 1, * } n , where components having value * are said free, otherwise they are fixed. The number of free components in a sub-hypercube h ∈ {0, 1, * } n is denoted by

rank(h) = |{i ∈ {1, . . . , n} | h i = * }|. A sub-hypercube h has 2 rank(h) ver- tices, denoted by c(h) = {x ∈ B n | ∀i ∈ {1, . . . , n}, h i ∈ B ⇒ x i = h i }. A sub-hypercube h ∈ {0, 1, * } n is smaller than a sub-hypercube h ∈ {0, 1, * } n whenever for each i ∈ {1, . . . , n} fixed in h (h i ∈ B), it is fixed to the value in h (h i = h i ). Then, remark that c(h ) ⊆ c(h).
A sub-hypercube h is closed by a BN f if the result of f applied to any of its vertices is one of its vertices: ∀x ∈ c(h), f (x) ∈ c(h); it is also known as a trap space. Given components K ⊆ {1, . . . , n}, h is K-closed by f whenever, for each component i ∈ K, either i is free in h, or f i applied on any vertices of h results in the fixed value h i . In other words, for all configurations in the K-closed sub-hypercube h, the next states of the components i ∈ K are in h:

∀x ∈ c(h), ∀i ∈ K, h i = * ⇒ f i (x) = h i . (1)
Let us consider the BN f defined in (B). The sub-hypercube h = 0 * * is closed by f , where c(0 * * ) = {000, 010, 001, 011}. This trap space indicates that component 1 can never get activated once deactivated. The sub-hypercube h = 0 * 1 is {1, 2}-closed by f as f 1 (001) = f 1 (011) = 0. Sub-hypercubes 1 * 0 and * * 0 are also {1, 2}-closed by f , where 1 * 0 is smaller than * * 0. In contrast, 1 * 1 is not {1, 2}-closed by f .

The Most Permissive update mode

Whenever a component changes state, e.g., increases from its minimal value 0 to its maximal value 1, the MP mode captures any behavior that could arise in the course of this change. For example, it may be that at some point, the component becomes high enough to activate one of its targets, whereas it remains not high enough to activate another one (because it has not reached its maximal value yet). These dynamic states can be captured by sub-hypercubes, where the changing components are free: they can be read both as 0 or 1.

Formally, the MP update mode can be defined as follows. Given a set of components K ⊆ {1, . . . , n}, let us denote by h x,K the smallest sub-hypercube of dimension n that contains x and that is K-closed by f . There is an MP transition from x to y whenever (a) y is a vertex of h x,K , and (b) the new state of all the components in K can be computed from h x,K :

∀x, y ∈ B n , x → MP y ⇐⇒ ∃K ⊆ {1, . . . , n} : y ∈ c(h x,K ) ∧ ∀i ∈ K, ∃z ∈ c(h x,K ) : y i = f i (z) .
h 001,{1,2} = * * 1 h 001,{1,2,3} = * * 1 h 011,{2,3} = 0 * 1 Fig. 2. Computation of sub-hypercubes with BN f (x) = (1, x1, (¬x1 ∧ x2) ∨ x3)
The abstraction properties and complexity results are detailed in [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF]. One can remark that the transition and reachability relations are identical: y is reachable from x if and only if x → MP y.

As this formalism serves as the cornerstone for the simulation algorithm, let us illustrate it by using the BN (A). Fig. 2 depicts the smallest sub-hypercubes with respect to a configuration x and a set K. Bold configurations belong to h x,K and those satisfying reachability condition (b) are boxed. The left side focuses on the initial configuration 001 and shows an iterative computation of

h 001,K from K = ∅ to K = {1, 2, 3}. As component 1 can flip in configuration 001, one obtains h 001,{1} = * 01. If component 1 is active, component 2 can flip, leading to h 001,{1,2} = * * 1. For all configurations z ∈ c h 001,{1,2} , component 3 stays active: f 3 (z) = 1. Hence, one obtains h 001,{1,2,3} = * * 1.
Here, some reachable configurations from x can be deduced. For instance, transitions 001 → MP 101 and 001 → MP 111 exist because reachability properties (a) and (b) are verified, However, 000 and 011 do not verify properties (a) and (b) respectively. Multiple sub-hypercubes may be required to capture all the MP transitions. This is illustrated in Fig. 2 on the right side with the configuration 011. The sub-hypercube h 011,{1,2,3} does not capture the transition to 001 because property (b) is not satisfied, while h 011,{2,3} does.

3 Simulation algorithm

Main principle

In its basic formal definition (2), the MP update mode defines the possible next configurations by considering all the subsets of components K ⊆ {1, . . . , n} and compute the associated sub-hypercube h x,K , that is the smallest K-closed subhypercube containing x. However, part of these subsets may be redundant. Let H be the components free in h x,K (3). By the minimality of h x,K , it means that for each component i ∈ H, there exists a configuration y ∈ c(h x,K ) such that f i (y) = x i . The set H can then be split in two: the components J for which their local function can be evaluated both to 0 and 1 from the configurations of h x,K (4), and the components L for which their local function is always evaluated to ¬x i from any of the configurations of h x,K :

H = i ∈ K | h x,K i = * , (3)
J = i ∈ H | ∃y, z ∈ c(h x,K ) : f i (y) = 0 ∧ f i (z) = 1 , (4)
L = i ∈ H | ∀y ∈ c(h x,K ) : f i (y) = x i . (5)
We have H = J ∪ L and J ∩ L = ∅. Then, remark that in (2), for a fixed x and K, there exists an MP transition x → MP y if and only if L ⊆ ∆(x, y) ⊆ H. Indeed, for each component i ∈ {1, . . . , n}, if y i = x i then necessarily i ∈ H. Moreover, if i ∈ L, the second condition of (2) imposes that y i = ¬x i . Whenever L = ∅, it results that for any strict subset K K, the transitions generated from the sub-hypercube h x,K form a (strict) subset of transitions generated from h x,K : the set of free components is a (strict) subset of H. Therefore, it is useless to explore any subset of K.

Whenever L = ∅, by the K-closeness property of h x,K , changing the state of a component i ∈ L would be irreversible while updating only components in K. Let us illustrate it by using the previous example without the ∨x 3 part:

f (x) = (1, x 1 , ¬x 1 ∧ x 2 )
. With the initial configuration x = 001, one can obtain the sub-hypercube h 001,{1,2,3} = * * * , L = {1} and J = {2, 3}. L = {1} means that component 1 cannot return to its initial state when it has begun to flip. Moreover, by the definition of MP transitions, all the components in L are modified. Therefore, whenever L is not empty, one should consider all the subhypercubes h x,K with (H \ L) ⊆ K H, to account for all these potential dependencies. Let us consider K = H \ L = {2, 3} in the previous example. One can obtain h 001,{2,3} = 00 * and L = {3}, which allows discovering a new transition: 001 → MP 000.

This approach leads us to compute a set of sub-hypercubes h x,K from subsets K of {1, . . . , n}, and where each of them can be characterized by a triplet (x, H, L) to which we refer to as a space. Each space (x, H, L) characterizes a set of MP transitions from x where the state of all the components in L is flipped, as well as any subset of components in H \ L: the transitions x → MP y with L ⊆ ∆(x, y) ⊆ H. In our algorithm, the set of transitions generated by a space is never enumerated explicitly, as it is exponential in |H \ L|. Instead, we count the number of transitions changing m components, from m = |L| (or 1 whenever L is empty) to |H|, i.e., the binomial coefficient |H\L| m-|L| . Overall, the sampling of the configuration following x can be summarized with the following steps:

1. compute a set of spaces S = {(x, H 1 , L 1 ), . . . , (x, H q , L q )}, corresponding to the sub-hypercubes {h x,K 1 , . . . , h x,K q }: starting with K 1 = {1, . . . , n}, we consider the sub-hypercubes closed along K 1 \ with each ⊆ L 1 , and recursively (thus whenever L 1 = ∅, q = 1). 2. for each space, count the number of transitions it can generate per number of changing components; 3. generate one random number to select the space (x, H, L) and the number m of components to flip; 4. randomly select m-|L| components in H \L, let us denote by C these chosen components; 5. flip the state of the components in L ∪ C.

The size of sets H and L is at most n. However, the number of spaces q can be exponential with n in the worst case.

Now, let us focus on the computation of h x,K , the smallest sub-hypercube containing x and K-closed by f . The main principle is to start from the subhypercube of rank 0 with x as the sole vertex, and iteratively free components to fulfill the K-closure property. To do so, we collect the set of components (among K) which can flip of state from at least one vertex of the sub-hypercube: for each i ∈ K, if there exists a vertex y of the sub-hypercube such that f i (y) = x i , then i must be free to verify the closeness property. This process is then repeated until the K-closeness property is verified, which in the worst case may require n iterations. It appears that the transitions generated by the asynchronous update mode match with the components marked as free in the first iteration only: the components i ∈ {1, . . . , n} such that f i (x) = x i . Thus, the additional transitions brought by the MP update mode are discovered in the later iterations only. This highlights the concept of permissive depth we introduce in this paper: the number of iterations in the computation of h x,K required to discover an MP transition. A depth of 1 corresponds to the asynchronous transitions, while a depth of n corresponds to the full MP dynamics. The simulation algorithm we propose in this paper allows controlling the depth of MP transitions, for instance by following a probabilistic distribution.

Algorithm

Listings 1 and 2 detail the steps for computing the reachable spaces and sampling the next configurations from them, as sketched in the previous section. In the description of the algorithms, we assumed fixed (1) a BN f of dimension n;

(2) a function depth to determine the depth threshold for computing the subhypercubes (it has to be an integer between 1 and n). For instance, it can be The spread function computes the K-closure of the sub-hypercube starting from the configuration x and stops after d iterations, d being a given depth. Determining whether the component i can change state requires determining the existence of a vertex z of the sub-hypercube such that f i (z) = x i . This is an The given algorithm makes the assumption that the BN f is locally monotone. In that case, given a sub-hypercube h, one can build in linear time from the influence graph G(f ) a configuration x min ∈ c(h) so that f i (x min ) = min{f i (x) | x ∈ c(h)}: the idea is that for each component j free in h, if j has a positive (monotone) influence on i (i.e., j +1 --→ i ∈ G(f )), then x min j = 0; if it has a negative influence on i, then x min j = 1. If j has no influence on i, its state in x min j can be arbitrary. A similar reasoning can be applied to compute a configuration x max ∈ c(h) so that f i (x max ) = max{f i (x) | x ∈ c(h)}. Then, one can decide in linear time whether i can change state in x in the scope of the sub-hypercube h: if x i is 1, then f i (x min ) must be 0, and if x i is 0, then f i (x max ) must be 1. Without this assumption, the can_flip function should be replaced with a call to a SAT solver.

H = H ∪ {i for i in K if can_flip(x, i, H, x[i])} 14 K = K \ H

Correctness, complexity, and parametrization

The sampling of the next configuration is driven by two parameters: a distribution over permissive depth, and a weight W ∈ R n ≥0 for transitions depending on the number of components they flip. These parameters can affect the generated dynamics and the complexity of the sampling.

General case: full MP dynamics. Let us assume that the depth function can always return n (either it is a constant function returning n, or the returned values follow a discrete distribution where the probability of drawing n is not 0), and that W ∈ R n >0 . Then, any MP transition has a non-zero probability to be sampled.

Given a space (x, H, L), as computed by reachable_spaces and characterized by a configuration x, a set H ⊆ {1, . . . , n} of free components, and a subset L ⊆ H of irreversible components, let us denote by tr(x, H, L) the set of candidate next configurations considered by sample_next_configuration:

tr(x, H, L) = {y ∈ B n | L ⊆ ∆(x, y) ⊆ H} . ( 6
)
We prove that the set of transitions the algorithm can generate from the computed spaces is equal to the full MP dynamics, except self-loops (Lemma 1), and that spaces generate disjoint sets of transitions (Lemma 2).

Lemma 1. Given a BN f of dimension n and one of its configuration x ∈ B n , and denoting by S the set of spaces returned by reachable_spaces(x) function,

(x,H,L)∈S tr(x, H, L) = {y | x → MP y, x = y} .
Lemma 2. Given a BN f of dimension n and one of its configuration x ∈ B n , and denoting by S the set of spaces returned by reachable_spaces(x), for any distinct pair of spaces (x, H, L), (x, H , L ) ∈ S, tr(x, H, L) ∩ tr(x, H , L ) = ∅.

Proofs are given in Appendix A.

Worst-case complexity. Assuming locally-monotone BNs, can_flip is performed in linear time (we assume the influence graph G(f ) is given); spread makes in the worst case n 2 call to can_flip, resulting in a cubic time in n; whereas irreversible is quadratic in n. Function reachable_spaces can then generate an exponential number of spaces. The sampling is then linear with respect to the number of spaces generated, and linear with respect to n. In the non-monotone case, can_flip is an NP-complete problem which currently can be solved in exponential time and space with SAT solvers.

Unitary depth: asynchronous and fully-asynchronous dynamics Let us consider the case whenever depth function always returns 1. The algorithm computes only one space (x, H, L = ∅) with H being the set of components i such that f i (x) = x i , and can generate any transition to y = x where ∆(x, y) ⊆ H. This corresponds exactly to the (general) asynchronous dynamics → g assuming W ∈ R n >0 . Moreover, as only one space is computed, the complexity drops to being linear in n, without any assumption on the local-monotony of f . Finally, whenever W 1 > 0 and for all m ∈ {2, . . . , n}, W m = 0, only transitions modifying one component are generated, matching exactly with the fullyasynchronous dynamics → a with equiprobable transitions, and with complexity similar to the previous restriction.

Sampling reachable attractors

The simulation of BNs is typically employed to assess the probability of reaching the different attractors. Because determining whether a configuration belongs to a cyclic attractor is a PSPACE-complete problem with (a)synchronous update modes [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF], most simulation algorithms are parametrized with a maximum number of steps to sample, without guarantee that an attractor has been reached. In the case of the MP update mode, the attractors turn out to be exactly the smallest closed sub-hypercubes (minimal trap spaces) of f , which can be computed at a much lower cost, albeit still relying on SAT solving [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF]. Thus, instead of fixing an arbitrary number of simulation steps, one can verify during the simulation whether the current configuration belongs to an attractor, and stop the sampling in that case. (sample_reachable_attractor of Listing 3).

In practice, the number of attractors reachable from a fixed initial configuration is usually small, and can be efficiently enumerated with the MP update mode, for instance using the mpbn tool [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF] 1 . In that case, the MP simulations can be employed to estimate the probability of reaching the different attractors depending on the depth and weight parameters. We then proceed as follows. Before simulating, we first compute the full set of attractors reachable from the initial configuration x. This set is then progressively refined during the simulation by removing attractors which are no longer included in h x,{1,...,n} . The simulation can then stop as soon as only one attractor can be reached, i.e., the current configuration is in the strong basin of an attractor (sample_reachable_attractor_bis of Listing 3).

These two schemes assume that the full MP dynamics is sampled, or that the attractors of the sampled dynamics match with the MP attractors as follows: each MP attractor is a superset of one and only one attractor of the sampled dynamics, and each attractor of the sampled dynamics is a subset of one and only one MP attractor. This is always the case for BNs having no cyclic attractors.

Evaluation

Using a prototype written in Python2 , we demonstrate the applicability of MP simulation on several Boolean models from the literature comprising about thirty components. After an illustration on toy examples, we study the effect of the parametrization on the assessment of the robustness of mutations impacting the propensities of reachable attractors in those models.

Our MP simulation algorithm has two parameters: the permissive depth and the weight of transitions W depending on the number of binary state changes. For the permissive depth, we will consider the constant n (full MP dynamics), constant 1 (asynchronous dynamics), and random sampling from a discrete exponentially decreasing distribution: depth d has probability 1/(2 d .M ) with M = n i=1 1/2 i the normalization factor. As depth n has a non-zero probability, this parametrization also enables the full MP dynamics, although largely prioritize transitions from low permissive depths. Regarding W , we will consider the uniform weight 1 n (random walk), and one-change only 10 n-1 .

Toy examples

Let us first consider the bi-stable example (A) from the configuration 000. This BN has two fixed points: 110 and 111. As explained in Sect. 1, the (a)synchronous dynamics from 000 predicts only one trajectory: 000 → g 100 → g 110. Thus, only one attractor is reachable (100% propensity). The MP dynamics uncovers another reachable attractor (111), as depicted in Fig. 3(left). Considering all the transitions, the MP simulation would conclude that the reachability of these two attractors are equiprobable: from the 4 initial transitions, 2 lead to the strong basin of 110 and two to the strong basin of 111. Whenever the depth is a random variable with an exponentially decreasing distribution, reaching 111 requires drawing a depth of 3 (probability of 1/7), leading to reaching 111 with probability 1/14, and 110 with probability 13/14. This indicates the sensitivity of the reachability of 111 to the permissive depth, and thus to the time scales of the underlying quantitative model.

Let us now consider the example (B) where the two attractors reachable with the asynchronous mode (Fig. 1 mode (Fig. 3(right)). In fully-asynchronous mode, they have an equal propensity, however in the (general) asynchronous mode, 001 propensity is twice the one of 110. Whereas the MP dynamics uncover additional trajectories and reachable configurations, in this particular example, the propensities of attractors is the same as with the general asynchronous case. An important aspect of MP dynamics is that it is transitive: any state reachable in a sequence of transitions is reachable in one transition. Thus, the size of the (reachable) strong basins of attractors can contribute significantly to their propensities. In the above two examples, the propensities of the attractors in MP uniform random walk is entirely determined by the size of their basins.

It should be also stressed that the parameters for determining the transition probabilities are usually arbitrary. Thus, analyzing the impact of parameters on the probability of reaching the different attractors by a random walk in the MP dynamics may only give an empirical insight on the Boolean dynamics, and cannot be formally transfered to an associated quantitative model.

Models from literature with different mutation conditions

From the literature, we selected BNs modeling cell fate decision processes: the reduced cell death receptor model [START_REF] Calzone | Mathematical modelling of cell-fate decision in response to death receptor engagement[END_REF] with 14 components; the tumor invasion model [START_REF] Cohen | Mathematical modelling of molecular pathways enabling tumour cell invasion and migration[END_REF] with 32 components; and the bladder model [START_REF] Remy | A modeling approach to explain mutually exclusive and co-occurring genetic alterations in bladder tumorigenesis[END_REF] with 35 components. These models have been designed with the fully-asynchronous update mode, and have been evaluated with respect to their ability to predict the changes of the attractors propensities subject to different mutation conditions (modeled by forcing some components to some state). For each model, we performed 10,000 simulations from the relevant initial configurations and for several mutation conditions, with different simulation parameters. With our prototype, the computational cost of the permissive depth is substantial, the simulations being between 3 times to 50 times slower than with depth 1 only, depending on the number of reachable spaces computed at each simulation step. However, no particular optimization has been implemented. Fig. 4 gives an example of the estimated propensities of attractors reachable from a unique initial configuration of the tumor invasion model with different mutation conditions and parameters. On the one hand, we can observe that the permissive depth can change drastically the proportions of reachable attractors: there are much more paths to the Metastasis attractor when considering potential time scale heterogeneity. The difference can be reduced by still considering the full MP dynamics, but giving exponentially decreasing probability to transitions with high permissive depth (computing transitions at depth at least D is done with probability n d=D (2 d • n i=1 2 -i ) -1 , with n the number of nodes in the model, 32 in this case). On the other hand, the qualitative effect of the p53 mutation is similar whenever analyzed only with the asynchronous simulations (permissive depth = 1), or in the most permissive cases: the propensity of the Metastasis attractor is larger in the mutant.

From a global analysis on the 3 models 3 , we observe that the absolute value of the propensities of the reachable attractors changes substantially when considering MP dynamics, in particular with depth fixed to n and uniform transition rates. Variable depth enables reducing the difference with the asynchronous dynamics while giving access to the full MP dynamics. Regarding the effect of mutations, it appears that, in most cases, the propensities of reachable attractors are affected in qualitatively the same direction, with various MP simulation parameters. This may indicate that these perturbations should be robust to heterogeneous time scales in the quantitative models captured by the BN abstraction. On the contrary, perturbations having a different qualitative effect on attractor propensities when considering permissive depth may be sensitive to the variety of time scales of the actual system, which could not be captured with usual (a)synchronous modes.

Discussion

The simulation of (asynchronous) BNs is a usual task in systems biology applications for assessing the effect of genetic perturbations on the propensities of the different cellular phenotypes. However, it is known that the asynchronous dynamics of BNs can preclude behaviors observed in quantitative systems, which may lead to substantial biases for the aforementioned analysis. On the contrary, the MP update mode ensures the completeness of the Boolean dynamics.

In this paper, we provide a first algorithm to sample trajectories from MP dynamics. The additional transitions predicted by MP come from iterative computations of possible state changes for each component, and where the first iteration corresponds to the asynchronous updates. By parameterizing the depth of this computation, we obtain a generalization of transition sampling from BNs with MP, (general) asynchronous, and fully-asynchronous dynamics. The main bottleneck of the provided algorithm is its potential exponential blow-up when considering all reachable spaces. Further work will address its efficient implementations and approximations.

As illustrated on different models from the literature, the simulation with the MP mode can largely affect predictions, both in terms of reachable attractors, and in terms of their propensities. Because of the transitive properties of MP dynamics (there is a direct transition to any reachable configuration), the attractors having a large (reachable) strong basin will dominate an equiprobable random walk of MP dynamics. Variable-depth simulation can then give more weight to transition with low permissive depth and sooth the effect of the size of the basins while still capturing the full MP dynamics.

It should be noted that, as with usual simulations of asynchronous BNs [START_REF] Mendes | Estimating attractor reachability in asynchronous logical models[END_REF][START_REF] Stoll | Continuous time boolean modeling for biological signaling: Application of gillespie algorithm[END_REF], the estimated attractor propensities are purely empirical and do not relate formally to the actual propensities in the modeled quantitative system. Due to its abstraction level, a BN is intrinsically non-deterministic. Assigning probabilities to transitions is a very strong assumption which, in this context, cannot be justified with any modeling or physical principle. At this level of abstraction, the behavior of the modeled system is not a Markov process and cannot be approximated by a Markov process (would the coefficient be fitted on data, for instance).

From a modeling perspective, our algorithm could be extended to have the permissive depth being different among components, enabling to fine-tune the time scale of their updates: the slower components should allow the more permissive depth. This would bring further control over the transitions added by the MP mode, while capturing trajectories precluded by the asynchronous mode.

A Proofs

A.1 Proof of Lemma 1 Consider any (x, H, L) ∈ S. Let h ∈ {0, 1, * } n be the sub-hypercube where dimensions in H are free, and otherwise fixed as in x: for each i ∈ {1, . . . , n}, if i ∈ H, then h i = * , otherwise, h i = x i . In such a configuration y, for each i ∈ H, there exists a configuration z ∈ c(h) such that f i (z) = y i : indeed, if y i = x i , then i / ∈ L. Thus, x → MP y by instantiating Eq. ( 2) with K = H. Conversely, let us assume that there exists K ⊆ {1, . . . , n} leading to x → MP y, and let H 0 be the set of free components in h x,K (note that ∆(x, y) ⊆ H 0 ). Let (x, H 1 , L 1 ) be the first space computed by our algorithm (i.e., from K = {1, . . . , n}). By construction, H 1 is the set of free components in h x,{1,...,n} , thus H 0 ⊆ H 1 . Two cases arise: either L 1 ⊆ ∆(x, y) and then the transition to y is generated from this first space, or there exists at least one component i ∈ L 1 whereas x i = y i . Recall that a component i can be in L 1 only if from each vertex z of the sub-hypercube f i (z) = x i . Thus K is necessarily a strict subset of H 1 , excluding at least the component i. Therefore, H 0 ⊆ K ⊆ H 1 \ {i ∈ L 1 | x i = y i } = H 2 , and by construction, (x, H 2 , L 2 ) ∈ S with L 2 ⊆ ∆(x, y).

A.2 Proof of Lemma 2

The set of spaces computed by reachable_spaces can be inductively characterized by a map S l for some l ∈ N from sets of components K to their associated subhypercube h x,K characterized by (H, L), as defined below (x being fixed it is omitted): Let us first consider the cases whenever L ⊆ H or L ⊆ H: by tr definition, tr(x, H, L) ∩ tr(x, H , L ) = ∅. Indeed, in the first case, remark that ∀y ∈ tr(x, H, L), L\H ⊆ ∆(x, y) while ∀y ∈ tr(x, H , L ), (L\H )∩∆(x, y ) = ∅, thus y = y . The second case is a symmetry.

We establish the following propositions:

- Consider any K → (H, L), K → (H , L ) ∈ S k , K = K , such that both L ⊆ H and L ⊆ H . By P4, there exists K → (H , L ) ∈ S k-1 with (K ∪K ) ⊆ K and such that there exists i ∈ L where i / ∈ K (note: it always work with K = {1, . . . , n}). If i ∈ K , then i ∈ L (by P3), thus L ⊆ H, a contradiction. Thus, i / ∈ K . By induction using P2 and P4, we obtain that K ⊆ K . By symmetry (apply the same reasoning by swapping K and K ), K = K .

Fig. 1 .

 1 Fig. 1. Transitions of the BN (B) from the configuration 111 with different modes

h

  001,∅ = 001 h 001,{1} = * 01 h 011,{1,2,3} = * * 1

Listing 1 . 2 # assumes f is locally monotone 3 if

 123 Computation of reachable spaces with the MP update mode 1 def can_flip(x: configuration, i: index, H: index set, v: bool):

8 9 10 #

 810 def spread(x: configuration, K: index set, d: depth): returns subset of K that can flip within the given depth 11

13

 13

16 17

 16 def irreversible(x: configuration, H: index set): 18 # returns subset of H that cannot flip back 19 return {i for i in H if not can_flip(x, i, H, 1-x[i])} 20 21 def reachable_spaces(x: configuration):22 d = depth() 23 S = [] # map of index set -> (index set, index set) (x, K, d) # H is subset of K 29 L = irreversible(x, H) if d > 1 else {} # L is subset of H 30for each M non-empty subset of L:31 J = K \ M 32 if J not in S and J not in Q: 33 Q.push(J) 34 S[K] = (H, L) 35 return Sa constant, or a sampling from a discrete distribution; (3) a weighting vector W ∈ R n ≥0 : W m is the weight of a transition modifying the state of m components simultaneously. A uniform random walk along MP transitions is thus obtained with W = 1 n and depth=n.

Listing 2 . 2 S 4 # compute apparent rate of transitions 5 R 7 if |L| > 0: 8 R

 224578 Sampling of the next configurations with the MP update mode 1 def sample_next_configuration(x: configuration) = reachable_spaces(x) 3 if |S| = 0: return x # fixed point = 0(len(S),n) # len(S) * n zero-filled matrix 6 for i = 1 to len(S): [i,|L|] = |L| * W[|L|] 9 for j = 1 to |H\ L|: 10 R[i,|L|+j] = binom(|H\L|,j) * W[|L|+j] 11 r = U[0, sum(R)[ # uniform sampling between 0 and sum(R) excluded 12 s,m = where(cumsum(R) > r) # s = space, m = nb of components to flip 13 H, L = S[s] 14 C = L ∪ random.sample(H\L, m -|L|) 15 y = copy(x) 16 y[C] = 1 -y[C] 17 return y instance of the classical Boolean satisfiability (SAT) problem, which, in general, is NP-complete.

Listing 3 .

 3 Algorithms for sampling a reachable attractor with the MP update mode 1 def sample_reachable_attractor(x: configuration):

4 x 6 if

 46 = sample_next_configuration(x) 5 every k iterations: # no need to verify at each step 1] # A contains only one element 10 11 def filter_reachable_attractors(A: sub-hypercube list, x: configuration): 12 H = spread(x, {1, . . . , n}, n) 13 return [a for a in A if a x/H] # a is smaller than the sub-hypercube 14 # formed by x and H 15 def sample_reachable_attractor_bis(x: configuration): 16 A = reachable_attractors(x) # list of attractors 17 while len(A) > 1: 18 x = sample_next_configuration(x) 19 every k iterations: # no need to verify at each step 20 A = filter_reachable_attractors(A, x) 21 return A[1] # A contains only one element

Fig. 3 .

 3 Fig. 3. MP transitions of BN (A) from configuration 000 (left) and of BN (B) from configuration 111 (right). Plain transitions have permissive depth 1, dashed depth 2, and dotted depth 3.

Fig.

  Fig. Estimated propensities of reachable attractors from a unique configuration of tumor invasion model in two mutation conditions (no mutation; p53 forced to 0), and different simulation parameters.

S 1 = 8 )

 18 {{1, . . . , n} → (H 0 , L 0 )} (7)S k+1 = S k ∪ {K \ M → (H , L ) | K → (H, L) ∈ S i , ∅ M ⊆ L} (Recall that ∀K → (H, L) ∈ S k , L ⊆ H ⊆ K, and tr(x, H, L) = {y ∈ B n | L ⊆ ∆(x, y) ⊆ H}.We prove that for any k ∈ N, ∀K → (H, L), K → (H , L ) ∈ S k , K = K , tr(x, H, L) ∩ tr(x, H , L ) = ∅ .

  P1) Any component i ∈ {1, . . . , n} such that there existsK → (H, L) with i ∈ L verifies f i (x) = ¬x i . Thus, ∀K → (H , L ), i ∈ K =⇒ i ∈ H . -(P2) By definition of S k , K → (H, L) ∈ S k =⇒ ∀i ∈ {1, . . . , n} \ K there exists K → (H , L ) ∈ S k-1 with i ∈ L . Thus by P1, f i (x) = ¬x i .-(P3) By sub-hypercube K-closeness definition and minimality,∀K → (H, L), K → (H , L ) ∈ S k , K ⊆ K =⇒ K ∩ L ⊆ L ⊆ H . -(P4) for any K = {1, . . . , n}, K → (H, L) ∈ S k \ S k-1 =⇒ ∃K →(H , L ) ∈ S k-1 with K K and K \ K ⊆ L (by P2).

https://github.com/bnediction/mpbn

https://github.com/bnediction/mpbn-sim

See https://doi.org/10.5281/zenodo.6725844 for supplementary information

Acknowledgments This work was supported by the French Agence Nationale pour la Recherche (ANR), with the project "BNeDiction" (ANR-20-CE45-0001).


Peter J Barry

Marcus A Mall

Antonio Álvarez

Carla Colombo

Karin M De Winter-De Groot

Isabelle Fajac

Kimberly A Mcbennett

Edward F Mckone

Bonnie W Ramsey

Sivagurunathan Sutharsan

P J Barry #

J L Taylor-Cousar

E Tullis

N Ahluwalia

L S Jun

S M Moskowitz

V Prieto-Centurion

S Tian

D Waltz

F Xuan

Y Zhang

S M Rowe #

D Polineni
email: dpolineni@kumc.edu

Vertex Pharmaceuticals

Triple Therapy for Cystic Fibrosis Phe508del-Gating and -Residual Function Genotypes

come   L'archive ouverte pluridisciplinaire

volume in 1 second (FEV 1 ) from baseline through week 8 in the elexacaftor-tezacaftor-ivacaftor group.

RESULTS-After the run-in period, 132 patients received elexacaftor-tezacaftor-ivacaftor and 126 received active control. Elexacaftor-tezacaftor-ivacaftor resulted in a percentage of predicted FEV 1 that was higher by 3.7 percentage points (95% confidence interval [CI], 2.8 to 4.6) relative to baseline and higher by 3.5 percentage points (95% CI, 2.2 to 4.7) relative to active control and a sweat chloride concentration that was lower by 22.3 mmol per liter (95% CI, 20.2 to 24.5) relative to baseline and lower by 23.1 mmol per liter (95% CI, 20.1 to 26.1) relative to active control (P<0.001 for all comparisons). The change from baseline in the Cystic Fibrosis Questionnaire-Revised respiratory domain score (range, 0 to 100, with higher scores indicating better quality of life) with elexacaftor-tezacaftor-ivacaftor was 10.3 points (95% CI, 8.0 to 12.7) and with active control was 1.6 points (95% CI, -0.8 to 4.1). The incidence of adverse events was similar in the two groups; adverse events led to treatment discontinuation in one patient (elevated aminotransferase level) in the elexacaftor-tezacaftor-ivacaftor group and in two patients (anxiety or depression and pulmonary exacerbation) in the active control group.

CONCLUSIONS-Elexacaftor-tezacaftor-ivacaftor was efficacious and safe in patients with

Phe508del-gating or Phe508del-residual function genotypes and conferred additional benefit relative to previous CFTR modulators. (Funded by Vertex Pharmaceuticals; VX18-445-104 ClinicalTrials.gov number, NCT04058353.) CYSTIC FIBROSIS IS A LIFE-SHORTENING autosomal recessive disease that affects more than 80,000 people worldwide. [START_REF] Bell | The future of cystic fibrosis care: a global perspective[END_REF][START_REF] Boeck | Progress in therapies for cystic fibrosis[END_REF][START_REF] Elborn | Cystic fibrosis[END_REF] In cystic fibrosis, deficiencies in the cystic fibrosis transmembrane conductance regulator (CFTR) protein, an epithelial anion channel, manifest as a complex multiorgan disease, including progressive respiratory impairment, exocrine pancreatic insufficiency, hepatobiliary disease, and abnormal sweat composition. [START_REF] Bell | The future of cystic fibrosis care: a global perspective[END_REF][START_REF] Elborn | Cystic fibrosis[END_REF] Measurement of sweat chloride concentration is a key diagnostic test for cystic fibrosis, and the concentration is used as an indicator of systemic CFTR function in clinical trials. [START_REF] Accurso | Sweat chloride as a biomarker of CFTR activity: proof of concept and ivacaftor clinical trial data[END_REF] Cystic fibrosis results from biallelic mutations in the CFTR gene. [START_REF] Bell | The future of cystic fibrosis care: a global perspective[END_REF] More than 1000 pathogenic CFTR mutations have been identified. [START_REF]Cystic Fibrosis Mutation Database (CFTR1)[END_REF][START_REF]The Clinical and Functional TRanslation of CFTR (CFTR2)[END_REF] Processing and trafficking mutations (e.g., Phe508del, the most common CFTR mutation) reduce the quantity of CFTR on the cell surface, [START_REF] Bell | The future of cystic fibrosis care: a global perspective[END_REF][START_REF]The Clinical and Functional TRanslation of CFTR (CFTR2)[END_REF][START_REF]Cystic Fibrosis Foundation patient registry: 2018 annual data report[END_REF] and channel-gating defects (e.g., Gly551Asp and other CFTR gating mutations) limit anion transport. CFTR mutations that result in lesser impairment of CFTR protein activity, collectively defined as residual function mutations, have also been identified. [START_REF] Bell | The future of cystic fibrosis care: a global perspective[END_REF][START_REF] Rowe | Tezacaftor-ivacaftor in residual-function heterozygotes with cystic fibrosis[END_REF] Most patients with gating or residual function CFTR mutations are heterozygous for the Phe508del mutation. [START_REF]The Clinical and Functional TRanslation of CFTR (CFTR2)[END_REF] Elucidation of the molecular consequences of CFTR mutations has supported the development of small-molecule modulators capable of restoring CFTR protein function. [START_REF] Van Goor | Rescue of CF airway epithelial cell function in vitro by a CFTR potentiator, VX-770[END_REF][START_REF] Ramsey | A CFTR potentiator in patients with cystic fibrosis and the G551D mutation[END_REF][START_REF] Van Goor | Nonclinical profile of the CFTR corrector VX-661[END_REF][START_REF] Taylor-Cousar | Tezacaftor-ivacaftor in patients with cystic fibrosis homozygous for Phe508del[END_REF][START_REF] Keating | VX-445-tezacaftor-ivacaftor in patients with cystic fibrosis and one or two Phe508del alleles[END_REF][START_REF] Middleton | Elexacaftor-tezacaftor-ivacaftor for cystic fibrosis with a single Phe508del allele[END_REF][START_REF] Mall | Cystic fibrosis: emergence of highly effective targeted therapeutics and potential clinical implications[END_REF] Ivacaftor, a CFTR potentiator, augments gating of mutant CFTR proteins. [START_REF] Van Goor | Rescue of CF airway epithelial cell function in vitro by a CFTR potentiator, VX-770[END_REF] In patients with gating mutations, ivacaftor improves lung function, nutritional status, and quality of life and decreases pulmonary exacerbations. [START_REF] Ramsey | A CFTR potentiator in patients with cystic fibrosis and the G551D mutation[END_REF][START_REF] Boeck | Efficacy and safety of ivacaftor in patients with cystic fibrosis and a non-G551D gating mutation[END_REF] Ivacaftor monotherapy is also efficacious and safe in patients with residual function mutations. [START_REF] Rowe | Tezacaftor-ivacaftor in residual-function heterozygotes with cystic fibrosis[END_REF] Tezacaftor is a CFTR corrector that ameliorates the defects in CFTR protein processing and cell-surface trafficking intrinsic to Phe508del. [START_REF] Van Goor | Nonclinical profile of the CFTR corrector VX-661[END_REF] Because Phe508del CFTR proteins also possess gating defects ,9 modulation requires both correction and potentiation. In patients heterozygous for Phe508del and specific residual function mutations, the combination of tezacaftor and ivacaftor improved lung function and sweat chloride concentrations as compared with ivacaftor alone. [START_REF] Rowe | Tezacaftor-ivacaftor in residual-function heterozygotes with cystic fibrosis[END_REF] The recently developed CFTR corrector elexacaftor has a mechanism of action that is complementary to that of tezacaftor. [START_REF] Keating | VX-445-tezacaftor-ivacaftor in patients with cystic fibrosis and one or two Phe508del alleles[END_REF] Pivotal studies showed that elexacaftor-tezacaftorivacaftor was efficacious and safe in patients with two Phe508del alleles and also in those with one Phe508del allele and an allele that makes no CFTR protein, which indicates that the presence of a single Phe508del allele is sufficient to confer responsiveness. [START_REF] Keating | VX-445-tezacaftor-ivacaftor in patients with cystic fibrosis and one or two Phe508del alleles[END_REF][START_REF] Middleton | Elexacaftor-tezacaftor-ivacaftor for cystic fibrosis with a single Phe508del allele[END_REF] These findings suggest that elexacaftor-tezacaftor-ivacaftor would provide additional clinical benefit in patients with Phe508del-gating and Phe508del-residual function genotypes by enhancing CFTR activity from the Phe508del allele. Here, we report results of a trial (VX18-445-104) designed to evaluate the magnitude of benefit of elexacaftor-tezacaftorivacaftor as compared with ivacaftor and tezacaftor-ivacaftor in patients 12 years of age or older with these genotypes.

METHODS

PATIENTS, TRIAL DESIGN, AND OVERSIGHT

This phase 3, multicenter, double-blind, parallel-group, randomized, active-controlled trial of elexacaftor-tezacaftor-ivacaftor enrolled patients 12 years of age or older with cystic fibrosis and Phe508del-gating or Phe508del-residual function genotypes. Complete inclusion and exclusion criteria as well as additional details on trial design, dosing, and statistical analysis are provided in the Supplementary Appendix, available with the full text of this article at NEJM.org; qualifying mutations are listed in Table S1.

To establish a reliable on-treatment baseline, patients entered a 4-week run-in period to receive either ivacaftor at a dose of 150 mg every 12 hours (ivacaftor comparator cohort; Phe508del-gating genotypes, including Phe508del-Arg117His) or tezacaftor at a dose of 100 mg once daily combined with ivacaftor (tezacaftor-ivacaftor comparator cohort; Phe508del-residual function genotypes) (Fig. S1). These genotype-defined comparator cohorts were based on approved indications for CFTR modulators in each country where the trial was conducted.

After the run-in period, patients who entered the 8-week treatment period were randomly assigned in a 1:1 ratio to receive either elexacaftor-tezacaftor-ivacaftor (elexacaftor, 200 mg once daily) or the regimen they received during the run-in period (ivacaftor or tezacaftorivacaftor). Stratification of randomization was determined during the run-in period on the basis of comparator cohort (ivacaftor vs. tezacaftor-ivacaftor), percentage of predicted forced expiratory volume in 1 second (FEV 1 ; <70 vs. ≥70), and sweat chloride concentration (<30 mmol per liter vs. ≥30 mmol per liter).

The trial was designed by Vertex Pharmaceuticals in collaboration with the authors. Each patient or the patient's legal guardian provided written informed consent, with assent obtained when age appropriate. Safety was monitored by an independent data monitoring committee. During the trial, the coronavirus disease 2019 pandemic led to the implementation of a global protocol addendum enabling patients to remain in the trial with in-home assessments and trialdrug provision. Data collection and analysis were performed by Vertex Pharmaceuticals in collaboration with the authors and the VX18-445-104 Study Group. The first two authors and last two authors wrote the first manuscript draft. All the authors had full access to the trial data after the final database lock, critically reviewed the manuscript, and approved it for submission. The investigators vouch for the accuracy and complete-ness of data generated at their respective sites, and the investigators and Vertex Pharmaceuticals vouch for the fidelity of the trial to the protocol, available at NEJM.org. Confidentiality agreements were in place between the sponsor and each investigative site during the trial.

END POINTS

The primary end point was the absolute change in the percentage of predicted FEV 1 from baseline through week 8 in the elexacaftor-tezacaftor-ivacaftor group. Key secondary end points, in hierarchical order, were the absolute change in sweat chloride concentration from baseline through week 8 in the elexacaftor-tezacaftor-ivacaftor group, the absolute change in the percentage of predicted FEV 1 from baseline through week 8 for elexacaftortezacaftor-ivacaftor as compared with active control (ivacaftor or tezacaftor-ivacaftor), and the absolute change in sweat chloride concentration from baseline through week 8 for elexacaftor-tezacaftor-ivacaftor as compared with active control. Other secondary end points were the absolute change in the score on the respiratory domain of the Cystic Fibrosis Questionnaire-Revised (CFQ-R; range, 0 to 100, with higher scores indicating a higher patient-reported quality of life with regard to respiratory symptoms) from baseline through week 8 in the elexacaftor-tezacaftor-ivacaftor group and for elexacaftor-tezacaftorivacaftor as compared with active control, as well as safety and the side-effect profile. The percentages of patients who reached sweat chloride concentrations below 60 mmol per liter and below 30 mmol per liter were assessed in a post hoc analysis.

STATISTICAL ANALYSIS

Efficacy analyses included all the patients who underwent randomization and received at least one dose during the treatment period. The absolute change from baseline in the percentage of predicted FEV 1 through week 8 was analyzed with the use of a mixed-effects model for repeated measures. The model included treatment group, visit, and treatmentgroup-by-visit interaction as fixed effects as well as continuous baseline percentage of predicted FEV 1 , continuous baseline sweat chloride concentration, and comparator cohort (ivacaftor vs. tezacaftor-ivacaftor) as covariates, with an unstructured covariance used for within-patient errors. The primary result that was obtained from the model was the estimated within-group change from baseline in the percentage of predicted FEV 1 through week 8 for elexacaftor-tezacaftor-ivacaftor. A similar mixed-effects model for repeated measures was applied to each of the key secondary end points. A hierarchical testing procedure was used to control the overall type I error rate at an alpha level of 0.05 for the primary and key secondary end points, the latter of which were prioritized in the testing hierarchy as a within-group analysis of sweat chloride concentration through week 8, a between-group analysis of the percentage of predicted FEV 1 through week 8, and a between-group analysis of sweat chloride concentration through week 8. For a hypothesis test to be considered statistically significant, the P value for that test and all the preceding tests in the hierarchy had to be below 0.05. Within-and between-group analyses of the other secondary end point of the absolute change in the CFQ-R respiratory domain score through week 8 were performed in a manner similar to the analyses of the primary and key secondary end points.

Subgroup analyses according to comparator cohort (ivacaftor [Phe508del-gating genotypes] vs. tezacaftor-ivacaftor [Phe508del-residual function genotypes]) were performed in a manner similar to the main analyses, including those for absolute changes in the percentage of predicted FEV 1 , in sweat chloride concentration, and in the CFQ-R respiratory domain score. Except for those involving the primary end point, these subgroup analyses were post hoc.

RESULTS

POPULATION

The trial was conducted at 96 sites in North America, Europe, and Australia, from August 28, 2019, to June 12, 2020. Overall, 271 patients entered the 4-week run-in period. After the run-in period, 258 patients (95 with Phe508del-gating genotypes and 163 with Phe508delresidual function genotypes) were randomly assigned to either the elexacaftor-tezacaftorivacaftor group (132 patients) or the active control group (126 patients) and received at least one dose of trial medication. Additional details on patient recruitment are provided in Figure S2. Treatment groups were well matched at baseline (Table 1 and Table S2). Demographic and clinical characteristics at baseline according to comparator cohort and treatment assignment (Table S3) and individual mutations on the second CFTR allele according to treatment assignment (Table S4) are provided in the Supplementary Appendix.

EFFICACY

Overall-The mean absolute change in the percentage of predicted FEV 1 from baseline (measured at the end of the run-in period) through week 8 with elexacaftor-tezacaftorivacaftor treatment (the primary end point) was 3.7 percentage points (95% confidence interval [CI], 2.8 to 4.6; P<0.001) (Table 2 and Fig. 1A). In contrast, the mean absolute change with active control (ivacaftor or tezacaftor-ivacaftor) was 0.2 percentage points (95% CI, -0.7 to 1.1), reflecting a between-group difference of 3.5 percentage points (95% CI, 2.2 to 4.7; P<0.001). The mean absolute change in sweat chloride concentration from baseline through week 8 with elexacaftor-tezacaftor-ivacaftor was -22.3 mmol per liter (95% CI, -24.5 to -20.2; P<0.001), as compared with a mean absolute change with active control of 0.7 mmol per liter (95% CI, -1.4 to 2.8), reflecting a between-group difference of -23.1 mmol per liter (95% CI, -26.1 to -20.1; P<0.001) (Table 2 and Fig. 1B). The mean absolute change in the CFQ-R respiratory domain score from baseline through week 8 with elexacaftor-tezacaftor-ivacaftor was 10.3 points (95% CI, 8.0 to 12.7) and with active control was 1.6 points (95% CI, -0.8 to 4.1), reflecting a between-group difference of 8.7 points (95% CI, 5.3 to 12.1) (Table 2 and Fig. 1C).

Subgroup Analyses-The results of a prespecified subgroup analysis of the primary end point (within-group absolute change from baseline in the percentage of predicted FEV 1 ) according to age at screening (<18 years vs. ≥18 years), sex, comparator cohort (Phe508delgating vs. Phe508del-residual function genotypes), percentage of predicted FEV 1 at baseline (<70 vs. ≥70), and geographic region (North America vs. Europe and Australia) were consistent with the results of the primary analysis (Fig. S3A). The same held for the results of a post hoc subgroup analysis of the between-group difference in the change from baseline in the percentage of predicted FEV 1 (Fig. S3B).

Additional post hoc analyses were performed to assess the role of genotype in treatment response. For patients with Phe508del-gating genotypes, including Phe508del-Arg117His, the mean absolute change in the percentage of predicted FEV 1 from baseline through week 8 with elexacaftor-tezacaftor-ivacaftor was 5.8 percentage points (95% CI, 4.2 to 7.4) and with ivacaftor control was 0.1 percentage points (95% CI, -1.6 to 1.7), for a difference of 5.8 percentage points (95% CI, 3.5 to 8.0); the mean absolute change in sweat chloride concentration with elexacaftor-tezacaftor-ivacaftor was -21.8 mmol per liter (95% CI, -25.7 to -17.8) and with ivacaftor was -1.8 mmol per liter (95% CI, -5.7 to 2.2), for a difference of -20.0 mmol per liter (95% CI, -25.4 to -14.6); and the mean absolute change in the CFQ-R respiratory domain score with elexacaftor-tezacaftor-ivacaftor was 10.2 points (95% CI, 6.6 to 13.8) and with ivacaftor was 1.3 points (95% CI, -2.5 to 5.2), for a difference of 8.9 points (95% CI, 3.8 to 14.0) (Table 3). The results of a post hoc analysis of data from patients with the Phe508del-Arg117His genotype were consistent with the results for the overall Phe508del-gating cohort (Table S5).

For patients with Phe508del-residual function genotypes, the mean absolute change in the percentage of predicted FEV 1 from baseline through week 8 with elexacaftor-tezacaftorivacaftor was 2.5 percentage points (95% CI, 1.4 to 3.5) and with tezacaftor-ivacaftor control was 0.5 percentage points (95% CI, -0.5 to 1.5), for a difference of 2.0 percentage points (95% CI, 0.5 to 3.4); the mean absolute change in sweat chloride concentration with elexacaftor-tezacaftor-ivacaftor was -23.1 mmol per liter (95% CI, -25.6 to -20.6) and with tezacaftor-ivacaftor was 1.7 mmol per liter (95% CI, -0.9 to 4.3), for a difference of -24.8 mmol per liter (95% CI, -28.4 to -21.2); and the mean absolute change in the CFQ-R respiratory domain score with elexacaftor-tezacaftor-ivacaftor was 10.4 points (95% CI, 7.2 to 13.7) and with tezacaftor-ivacaftor was 1.9 points (95% CI, -1.4 to 5.1), for a difference of 8.5 points (95% CI, 4.0 to 13.1) (Table 3).

After 8 weeks of treatment with elexacaftor-tezacaftor-ivacaftor, the mean (±SD) sweat chloride concentrations were 32.7±23.5 mmol per liter for patients with Phe508del-gating genotypes and 39.9±19.3 mmol per liter for those with Phe508del-residual function genotypes, as compared with 52.0±21.9 mmol per liter with ivacaftor and 63.4±27.3 mmol per liter with tezacaftor-ivacaftor for patients in those active control groups. An analysis of individual sweat chloride concentrations through week 8 for patients who received elexacaftor-tezacaftor-ivacaftor showed that 83.3% had concentrations below 60 mmol per liter and 50.0% had concentrations below 30 mmol per liter; for patients who received active control, the corresponding percentages were 55.5% and 17.6% (Fig. S4 and Table S6). For patients with Phe508del-gating genotypes, 65% who received elexacaftor-tezacaftorivacaftor were below 30 mmol per liter through week 8, as compared with 16% of those who received ivacaftor. For patients with Phe508del-residual function genotypes, 42% who received elexacaftor-tezacaftor-ivacaftor were below 30 mmol per liter through week 8, as compared with 19% of those who received tezacaftor-ivacaftor.

SAFETY

Overall, 66.7% of the patients in the elexacaftor-tezacaftor-ivacaftor group and 65.9% of those in the active control group (ivacaftor or tezacaftor-ivacaftor) had one or more adverse events, which for most patients were mild or moderate in severity and resolved during the trial (Table 4). Serious adverse events were reported in 5 patients (3.8%) in the elexacaftor-tezacaftor-ivacaftor group and 11 patients (8.7%) in the active control group, with the difference attributable to a higher incidence of pulmonary exacerbation in the active control group. One patient in the elexacaftor-tezacaftor-ivacaftor group discontinued treatment owing to an adverse event (elevated aminotransferase level), and 2 patients in the active control group discontinued treatment owing to an adverse event (anxiety or depression in 1 patient and pulmonary exacerbation in 1 patient).

On the basis of previous experience with elexacaftor-tezacaftor-ivacaftor, including the phase 3 trials, [START_REF] Middleton | Elexacaftor-tezacaftor-ivacaftor for cystic fibrosis with a single Phe508del allele[END_REF][START_REF] Heijerman | Efficacy and safety of the elexacaftor plus tezacaftor plus ivacaftor combination regimen in people with cystic fibrosis homozygous for the F508del mutation: a double-blind, randomised, phase 3 trial[END_REF] data regarding aminotransferase levels, rash, creatine kinase level, and blood pressure were reviewed. Elevated levels of alanine aminotransferase or aspartate aminotransferase that were greater than three times, greater than five times, and greater than eight times the upper limit of the normal range occurred in 4 of 125 patients for whom data were available (3.2%), 1 patient (0.8%), and 1 patient (0.8%), respectively, in the elexacaftor-tezacaftor-ivacaftor group and in 2 of 123 patients (1.6%), 1 patient (0.8%), and no patients, respectively, in the active control group (Table S7). Eight patients (6.1%) in the elexacaftor-tezacaftor-ivacaftor group and 1 patient (0.8%) in the active control group had adverse events involving elevated aminotransferase levels (Table 4). No patient had a serious adverse event involving elevated aminotransferase levels. Rash was observed in 4 patients (3.0%) in the elexacaftor-tezacaftor-ivacaftor group and in 5 patients (4.0%) in the active control group (Table S8). All cases of rash were mild or moderate in severity. Increased blood creatine kinase levels were reported in 2 patients (1.5%) in the elexacaftor-tezacaftorivacaftor group and in no patients in the active control group. Baseline mean systolic and diastolic blood pressures increased by 3.0 mm Hg and 2.5 mm Hg, respectively, in the elexacaftor-tezacaftor-ivacaftor group and by 0.5 mm Hg and 0.3 mm Hg, respectively, in the active control group at week 8 (Table S9). There were no adverse events involving hypertension in the elexacaftor-tezacaftor-ivacaftor group or the active control group. There were no notable safety findings in other clinical or laboratory assessments.

DISCUSSION

Here, we report results of an 8-week trial of elexacaftor-tezacaftor-ivacaftor in patients with cystic fibrosis and either Phe508del-gating or Phe508del-residual function genotypes.

Elexacaftor-tezacaftor-ivacaftor treatment improved lung function and sweat chloride concentration relative to an active control (ivacaftor or tezacaftor-ivacaftor). Most patients had adverse events that were mild or moderate in severity and that were consistent with those observed in previous studies. [START_REF] Middleton | Elexacaftor-tezacaftor-ivacaftor for cystic fibrosis with a single Phe508del allele[END_REF][START_REF] Heijerman | Efficacy and safety of the elexacaftor plus tezacaftor plus ivacaftor combination regimen in people with cystic fibrosis homozygous for the F508del mutation: a double-blind, randomised, phase 3 trial[END_REF] These results confirm that by enhancing activity from the Phe508del allele, elexacaftor-tezacaftor-ivacaftor can provide additional benefit to patients with a single Phe508del allele plus a gating or residual function allele that is responsive to previous CFTR modulator regimens.

For patients with CFTR gating mutations, the current benchmark for effective treatment is ivacaftor monotherapy, which has been shown to improve the percentage of predicted FEV 1 by 10.6 percentage points in patients with the Gly551Asp gating mutation relative to placebo. [START_REF] Ramsey | A CFTR potentiator in patients with cystic fibrosis and the G551D mutation[END_REF] In previous studies, tezacaftor-ivacaftor improved the percentage of predicted FEV 1 by 6.8 percentage points in patients with Phe508del-residual function genotypes relative to placebo. [START_REF] Rowe | Tezacaftor-ivacaftor in residual-function heterozygotes with cystic fibrosis[END_REF] Despite the clinical heterogeneity of these genotype groups, elexacaftor-tezacaftor-ivacaftor treatment resulted in additional increases in the percentage of predicted FEV 1 as compared with either ivacaftor or tezacaftor-ivacaftor treatment; this is largely attributable to enhanced function of CFTR protein arising from the Phe508del allele.14 In a post hoc subgroup analysis of the percentage of predicted FEV 1 , the treatment differences between elexacaftor-tezacaftor-ivacaftor and active control were 5.8 percentage points (95% CI, 3.5 to 8.0) for patients with Phe508del-gating genotypes and 2.0 percentage points (95% CI, 0.5 to 3.4) for patients with Phe508del-residual function genotypes. Relative to patients with Phe508del-gating genotypes, smaller changes in lung function in response to CFTR modulation among patients with Phe508del-residual function genotypes have been observed previously and probably reflect differences in the progression of lung disease in patients with Phe508del-residual function genotypes, who are generally older than those in other genotype groups. [START_REF] Rowe | Tezacaftor-ivacaftor in residual-function heterozygotes with cystic fibrosis[END_REF] (In the current trial, such patients had a mean age of 40.8 years, as compared with 32.2 years for patients with Phe508del-gating genotypes.)

Elexacaftor-tezacaftor-ivacaftor treatment also resulted in changes in CFQ-R respiratory domain scores, relative to ivacaftor or tezacaftor-ivacaftor, that exceeded the established minimal clinically important difference (4 points), a finding that shows further abatement of respiratory symptoms in these patients. [START_REF] Rowe | Tezacaftor-ivacaftor in residual-function heterozygotes with cystic fibrosis[END_REF][START_REF] Ramsey | A CFTR potentiator in patients with cystic fibrosis and the G551D mutation[END_REF][START_REF] Boeck | Efficacy and safety of ivacaftor in patients with cystic fibrosis and a non-G551D gating mutation[END_REF][START_REF] Quittner | Determination of the minimal clinically important difference scores for the Cystic Fibrosis Questionnaire-Revised respiratory symptom scale in two populations of patients with cystic fibrosis and chronic Pseudomonas aeruginosa airway infection[END_REF] Sweat chloride concentration may differentiate the effectiveness of CFTR modulator regimens at a population level. At baseline, the mean sweat chloride concentration in patients was approximately 60 mmol per liter, the threshold for definitive diagnosis of cystic fibrosis. [START_REF] Farrell | Diagnosis of cystic fibrosis: consensus guidelines from the Cystic Fibrosis Foundation[END_REF] After 8 weeks of elexacaftor-tezacaftor-ivacaftor treatment, 50.0% of the patients had sweat chloride concentrations below 30 mmol per liter, a level that matches those generally seen in the population of asymptomatic carriers with a single mutant CFTR allele and a level below which a diagnosis of cystic fibrosis is unlikely, [START_REF] Farrell | Diagnosis of cystic fibrosis: consensus guidelines from the Cystic Fibrosis Foundation[END_REF][START_REF] Wilschanski | Mutations in the cystic fibrosis transmembrane regulator gene and in vivo transepithelial potentials[END_REF] whereas only 17.6% of the patients who received active control had sweat chloride concentrations below 30 mmol per liter. This result reflects improved CFTR function with elexacaftor-tezacaftor-ivacaftor treatment.

This phase 3 trial showed the efficacy of elexacaftor-tezacaftor-ivacaftor therapy in patients with Phe508del-gating and Phe508del-residual function genotypes, with clinical benefit exceeding previous CFTR modulators. No new safety findings were noted as compared with previous studies involving patients with cystic fibrosis. [START_REF] Middleton | Elexacaftor-tezacaftor-ivacaftor for cystic fibrosis with a single Phe508del allele[END_REF][START_REF] Heijerman | Efficacy and safety of the elexacaftor plus tezacaftor plus ivacaftor combination regimen in people with cystic fibrosis homozygous for the F508del mutation: a double-blind, randomised, phase 3 trial[END_REF] Sweat chloride concentrations after elexacaftor-tezacaftor-ivacaftor treatment reached levels at or near those found in asymptomatic carriers with a single CFTR mutation. These findings confirm the efficacy and safety of elexacaftor-tezacaftor-ivacaftor in patients with at least one Phe508del allele. to a 100-point range, with higher scores indicating a higher patient-reported quality of life with regard to respiratory symptoms; the minimal clinically important difference (MCID) is 4 points and is indicated in the plot by the straight gray line. In Panels A through C, data are least-squares means, and the I bars indicate standard errors; the dashed line at 0 cor responds to no change from baseline. The sample size shown under each x axis is the number of patients at that time point with data that could be evaluated. ELX-TEZ-IVA denotes elexacaftor-tezacaftor-ivacaftor, IVA ivacaftor, and TEZ-IVA tezacaftor-ivacaftor.

‖

Scores on the respiratory domain of the Cystic Fibrosis Questionnaire-Revised (CFQ-R) range from 0 to 100, with higher scores indicating a higher patient-reported quality of life with regard to respiratory symptoms.

**

The body-mass index is the weight in kilograms divided by the square of the height in meters.

Figure 1 .

 1 Figure 1. Efficacy End Points. Panel A shows the absolute change from baseline at each visit in the percentage of predicted forced expiratory volume in 1 second (FEV 1 ) on the basis of a mixed-effects model for repeated measures. Panel B shows the absolute change from baseline at each visit in the sweat chloride concentration on the basis of a mixed-effects model for repeated measures. Panel C shows the absolute change from baseline at each visit in the score on the respiratory domain of the Cystic Fibrosis Questionnaire-Revised (CFQ-R) on the basis of a mixed-effects model for repeated measures. Respiratory domain scores are normalized

  multiple events within a category was counted only once in that category. † Infective pulmonary exacerbation of cystic fibrosis was the only serious adverse event that occurred in two or more patients in either treatment group. ‡ Shown are adverse events that occurred in at least 5% of the patients in either group. Adverse events are listed according to the preferred term in the Medical Dictionary for Regulatory Activities (MedDRA), version 23.0. N Engl J Med. Author manuscript; available in PMC 2022 April 05.

Table 1 .

 1 Demographic and Clinical Characteristics of the Patients at Baseline.

*

Characteristic

Elexacaftor-Tezacaftor-Ivacaftor (N = 132) Active Control: Ivacaftor or Tezacaftor-Ivacaftor (N = 126)

  Shown are the demographic and clinical characteristics of the full analysis set, which was defined as all randomly assigned patients who received at least one dose of trial drug in the treatment period. Baseline was defined as the most recent nonmissing measurement before the first dose of trial drug in the treatment period (after the 4-week run-in period). Percentages may not total 100 because of rounding. FEV 1 denotes forced expiratory volume in 1 second. † Age distribution was calculated on the basis of age at the time of screening. ‡ Ethnic group and race were reported by the patient. § There were no patients of Asian race or Native Hawaiian or other Pacific Islander ethnic group. If a patient was reported to have multiple races or ethnic groups, then the patient was counted for each race or ethnic group reported. ¶ Patients with Phe508del-gating genotypes who were assigned to the active control group received ivacaftor, and patients with Phe508del-residual function genotypes who were assigned to the active control group received tezacaftor-ivacaftor.

	Elexacaftor-Tezacaftor-Ivacaftor (N = 132) Active Control: Ivacaftor or Tezacaftor-Ivacaftor (N = 126)	2 (1.5) 2 (1.6)	70 (53.0) 63 (50.0)	53 (40.2) 52 (41.3)	7 (5.3) 9 (7.1)	65 (49.2) 65 (51.6) 59.5±27.0 56.4±25.5	67 (50.8) 61 (48.4) 76.5±16.6 77.3±15.8	24.07±4.72 24.05±4.71	37.7±14.7 37.6±14.3		15 (11.4) 9 (7.1)	117 (88.6) 117 (92.9)		5 (3.8) 4 (3.2)	117 (88.6) 114 (90.5)	10 (7.6) 8 (6.3)		122 (92.4) 111 (88.1)	0 2 (1.6)	0 1 (0.8)	1 (0.8) 4 (3.2)	9 (6.8) 9 (7.1)		49 (37.1) 48 (38.1)	70 (53.0) 64 (50.8)	13 (9.8) 14 (11.1)		50 (37.9) 45 (35.7)	82 (62.1) 81 (64.3)		67.1±15.7 68.1±16.4
	Characteristic	<40%	40 to <70%	70 to ≤90%	Sex -no. (%) >90%	Male Sweat chloride concentration -mmol/liter	Female CFQ-R respiratory domain score ‖	Age Body-mass index **	Mean -yr	Distribution -no. (%) † * Plus-minus values are means ±SD.	12 to <18 yr	≥18 yr	Hispanic or Latino ethnic group -no. (%) ‡	Hispanic or Latino	Not Hispanic or Latino	Not collected per local regulations	Race or ethnic group -no. (%) ‡ §	White	Black	American Indian or Alaska Native	Other	Not collected per local regulations	Geographic region -no. (%)	North America	Europe	Australia	Genotype -no. (%) ¶	Phe508del-gating	Phe508del-residual function	Percentage of predicted FEV 1	Mean	Distribution -no. (%)

N Engl J Med. Author manuscript; available in PMC 2022 April 05.

N Engl J Med. Author manuscript; available in PMC 2022 April 05.

from baseline through week 8 for elexacaftor-tezacaftor-ivacaftor as compared with active control, and the absolute change in sweat chloride concentration from baseline through week 8 for elexacaftor-tezacaftor-ivacaftor as compared with active control. N Engl J Med. Author manuscript; available in PMC

April 05.

Acknowledgments

Supported by Vertex Pharmaceuticals. The National Institute for Health Research (NIHR) provided support to the NIHR Manchester Clinical Research Facility. The National Institutes of Health provided support to the University of Alabama (P30DK072482 and UL1TR003096) and the University of Kansas Medical Center Research Institute (P20GM130423).

Disclosure forms provided by the authors are available at NEJM.org.

A data sharing statement provided by the authors is available with the full text of this article at NEJM.org.

We thank the patients and their families for participating in this trial and all site investigators and coordinators for their contributions to the trial; the members of the Cystic Fibrosis Foundation Therapeutics Development Network and the European Cystic Fibrosis Society Clinical Trials Network for their support of the trial sites; Emily Poulin, Ph.D., an employee of Vertex Pharmaceuticals, who may own stock or stock options in the company, for providing editorial coordination and support; and Nathan Blow, Ph.D., of Vertex Pharmaceuticals, who may own stock or stock options in the company, and Christopher Edwards, Ph.D., C.M.P.P., of ArticulateScience, for providing editorial assistance under the guidance of the authors and with support from Vertex Pharmaceuticals.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.

Appendix

The authors' full names and academic degrees are as follows: Peter J. Barry, M.D., Marcus A. Mall The key secondary end points, in hierarchical order, were the absolute change in sweat chloride concentration from baseline through week 8 in the elexacaftor-tezacaftor-ivacaftor group, the absolute change in the percentage of predicted FEV

Author Manuscript

Barry et al.

Table 3.

Subgroup Analysis of Efficacy Results