An integrated modelling approach for targeted degradation: insights on optimization, data requirements and PKPD predictions from semi- or fully-mechanistic models and exact steady state solutions

The value of an integrated mathematical modelling approach for protein degraders which combines the benefits of traditional turnover models and fully mechanistic models is presented. Firstly, we show how exact solutions of the mechanistic models of monovalent and bivalent degraders can provide insight on the role of each system parameter in driving the pharmacological response. We show how on/off binding rates and degradation rates are related to potency and maximal effect of monovalent degraders, and how such relationship can be used to suggest a compound optimization strategy. Even convoluted exact steady state solutions for bivalent degraders provide insight on the type of observations required to ensure the predictive capacity of a mechanistic approach. Specifically for PROTACs, the structure of the exact steady state solution suggests that the total remaining target at steady state, which is easily accessible experimentally, is insufficient to reconstruct the state of the whole system at equilibrium and observations on different species (such as binary/ternary complexes) are necessary. Secondly, global sensitivity analysis of fully mechanistic models for PROTACs suggests that both target and ligase baselines (actually, their ratio) are the major sources of variability in the response of non-cooperative systems, which speaks to the importance of characterizing their distribution in the target patient population. Finally, we propose a pragmatic modelling approach which incorporates the insights generated with fully mechanistic models into simpler turnover models to improve their predictive ability, hence enabling acceleration of drug discovery programs and increased probability of success in the clinic.


Introduction
Mechanistic modelling has proven extremely valuable not only in enhancing the understanding of both traditional and new modalities mechanism of action [1][2][3][4], but also in supporting and guiding compound optimization by shedding light on Structure-Activity Relationship (SAR). Especially for emerging new modalities where little is known about the mechanism or where technology lags behind science in generating reliable data on biological or pharmacological quantities of interest, inexpensive high-throughput model simulations can to some extent ''bridge the gap'' between science and technology by predicting the response of a biological system in scenarios of interest and, more importantly, by identifying which mechanistic parameters are key drivers of the response. Such quantitative understanding can ultimately provide a robust rationale to identify which missing data would be most informative in building an understanding of the pharmacology, and hence which technologies need prioritization for data generation.
On the other hand, though, the use of mechanistic models to explain available data can be impractical and potentially even misleading without a thorough understanding of the amount and type of data that is required to uniquely identify the model parameters. While modelling software is designed to always return parameter estimates, if the model structure is excessively articulated (overparametrized) for the data, or, conversely, if the amount and/or type of data is insufficient or too simplistic to inform the building blocks of a mechanistic model, multiple parameter estimates can provide an equally optimal data fit. However, the uncertainty on those estimates is typically high, or their reliability low, i.e. the output values are highly unlikely to be truly representative of the biological, chemical or pharmacological quantity they encode (e.g., protein baselines, endogenous turnover rates, on/off binding rates, drug-induced degradation rates, ...). In such case the risk is over-interpreting the parameter values, potentially leading to wrong decisions or scientific conclusions at different levels in a drug discovery cascade, from compound optimization to assessment of therapeutic potential.
Although much has been published on identifiability analysis of differential equations to address this issue [5,6], the minimal required data on target kinetics can be challenging to generate, and even when this is not the case its generation requires long timelines and significant resources. Actually, in many cases kinetics data provide a level of detail that goes beyond the necessary and sufficient needs for pragmatically understanding the underlying mechanism of action (MoA). On the contrary, the steady state (i.e., dynamical equilibrium) of a biological or pharmacological system can be enough to provide insights on the MoA and hence to build a robust and reliable predictive model from higher throughput and more easily accessible data [7][8][9][10]. Moreover, the steady state can often be mechanistically described by simpler (non-linear) algebraic equations directly derived from a parent set of differential equations, for which an exact solution might even be available (depending on the mathematical structure of the system) [7,8,11,12]. Although the resulting mathematical formulae may still retain some level of complexity (and they typically do), they can facilitate the identification of independent surrogate parameters, i.e. compositions of the original parameters into sums, ratios or other functional forms, which can lead to a reduction of the model parametrization and, ultimately, to improved model reliability. Such approach to identifiability can further be combined with global sensitivity analysis [13] to identify which of the original mechanistic parameters dominate the system response variability, and hence whose experimental measurements would be most impactful.
Whenever an exact steady state solution cannot be calculated or its complexity is still prohibitive, semi-mechanistic models such as indirect-response (turnover) models [14] offer an alternative option. Such models are easier to use in practice due to the smaller number of parameters that require fitting, however while they can successfully explain the data in specific studies or experimental settings, they lack information on binding kinetics or baseline levels, which makes them prone to failure at predicting the response for different compounds or cell lines [15,16].
In this manuscript we propose an integrated modelling approach which leverages the benefits of each type of model to mitigate the limitations of the others. Firstly we show how exact mechanistic steady state solutions of the MoA of binary-and ternary-complex degraders can (i) provide noise-free insight on the role of each system parameter in driving the pharmacological response regardless of its specific value, (ii) suggest which mechanistic knowledge can be confidently extracted from single time point data and (iii) which additional data needs to be collected to enhance the mechanistic understanding of the system, and (iv) ultimately, inform compound optimization, data generation and resource prioritization.
Secondly, we show how global sensitivity analysis of fully mechanistic models can help to identify the key drivers of the response.
Finally, we propose a pragmatic modelling approach which incorporates the insights generated with fully mechanistic models into simpler turnover models to improve their predictive ability, hence enabling acceleration of drug discovery programs and increased probability of success in the clinic.

Exact solution of bilinear systems
The life cycle of biological entities and their interaction with chemical compounds can be described via non-linear systems of ordinary differential equations (ODEs), where the type of non-linearity is often bilinear or at most quadratic [17]. Exact steady-state solutions of such systems are challenging to compute and rarely available in closed form, due to potentially many chemical species involved and corresponding interactions, resulting in a large number of non-linear equations. Nevertheless, in some cases implicit relationships can be obtained and solved numerically. Not only, they can inform on the system identifiability, i.e. which parameters (or surrogates thereof) can be uniquely estimated from the data. A mathematical method to obtain an implicit, exact steady state solution to chemical reaction networks with bilinear rate laws is described in [18] and is summarized in Appendix A. Briefly, such method applies thoughtful algebraic manipulations to leverage the linear component of the system while segregating the non-linear part to its core, which can then be solved numerically.

Monovalent degraders
In this paper we refer to compounds that induce degradation of a protein of interest (PoI) by binding solely to it as Monovalent Degraders (MDs). It is very likely that other components are required for degradation of the protein (e.g. the proteasome), nevertheless no assumption is made here on the mechanism of degradation. The only assumption is that the compound needs to bind solely to the PoI to induce its degradation with no requirement to bind any other component of the system (differently from bivalent degraders).

A mechanistic model for monovalent degraders
We assume that under endogenous conditions (i.e. in absence of compound) the target protein synthesis and degradation are zero-and first-order processes, respectively, with corresponding rates k syn , k deg . When compound is added, binding kinetics governed by on/off rates k T on , k T off leads to the formation of a binary complex, which induces degradation at a first order rate k MD . Upon degradation the compound is released and recycled (Fig. 1). Such mechanism can be mathematically described by the following system of ODEs: with initial conditions Tð0Þ ¼ T 0 ¼ k syn =k deg , Dð0Þ ¼ D 0 , TDð0Þ ¼ 0, where T, D, and TD represent target, drug, and binary complex time-varying concentrations, respectively. Being the MD recycled, the total MD concentration D 0 ¼ D þ TD is preserved. As a result, the last (or second) equation is redundant as the binary complex (or MD) concentration at any time can be obtained via the linear conservation law as Since we are interested in the steady state, we set the derivatives in (1) to 0 to obtain a steady state model: Note that, while in system (1) T, D and TD are timedependent variables, they are constant in (3) by definition of steady state (the same notation has been used for the sake of simplicity). Because system (3) is non-linear, calculating an exact solution is not straightforward. Nevertheless, this type of non-linearity (bilinear) falls within the category of tractable systems which can be solved with the mathematical method developed in [18].

PROTACs
Proteolysis targeting chimeras (PROTACs) are a novel drug modality that fosters degradation catalysis by coopting an E3 ligase (e.g. Cereblon or VHL) to tag the targeted PoI for turnover by the proteasome. Such mechanism of action heavily relies on the formation of a ternary complex with the PoI and E3 ligase, which is known to be liable to auto-inhibition, i.e., impairment of ternary complex formation at high PROTAC concentrations due to an excess of PoI-or E3 ligase-PROTAC binary complex [19,20]. In other words, while in a two-body binding system the amount of binary complex increases monotonically with the binder concentration up to ligand saturation following a sigmoidal relationship, ternary complex formation can decrease as the PROTAC concentration increases beyond a critical threshold and is hence typically described by a double sigmoidal (bell-shaped) function ( Fig. 2, left). From a pharmacological perspective, such binding dynamics in a biological setting can result into the so-called hook effect [19], i.e., a loss of degradation at high PROTAC concentrations, with the implication that maximal effect can only be achieved within a certain concentration window, whose center and breadth are highly specific to each molecule (Fig. 2, right). It is worth noting that not all PROTACs display auto-inhibition in practice. Efforts to reduce the hook effect whenever present have been summarized by Cecchini et al. [21], nevertheless it is fair to say that additional work is required to fully understand how to pragmatically minimize this phenomenon.

PROTACs mechanistic modelling
Ternary complex formation (and subsequent degradation) can happen through two different pathways: from a PRO-TAC-target (PT) or PROTAC-ligase (PL) binary complex, where the extent of the contribution of each pathway is dictated by binding affinities, PROTAC concentration, and target and E3 ligase baselines (Fig. 3). Once the ternary complex is formed the PoI is degraded while PROTAC and E3 ligase are released and recycled back into the system. It is important to keep in mind that PROTACs are not themselves degraders, rather degradation catalysts. Therefore, the apparent ''PROTAC degradation rate'' (k PRO ) is in fact a surrogate, composite parameter that synthesizes (i) ubiquitin transfer rate (which depends on the stereochemistry), (ii) ubiquitination rate and (iii) proteasomal degradation rate as a whole. While the stereochemistry can be optimized -at least in principle -to facilitate the ubiquitin transfer, the latter two parameters are endogenous to the biological system and can dictate the overall degradation kinetics (i.e. can be rate-limiting). This underscores the importance of understanding biological differences across cell lines, which cannot be conceived separately from the compound's kinetics.
In vitro data for different targets suggests that three degradation scenarios can occur to the PoI bound into a PT binary complex: 1. Endogenous degradation (flat line around the baseline) 2. Low-moderate degradation (more efficient than endogenous degradation, less efficient than ternary complex degradation), for example when the target warhead is a MD itself (shallow degradation curve) 3. Stabilization, i.e., the compound inhibits the endogenous degradation machinery (sigmoidal curve above baseline).
Such observation justified the introduction of a binary complex degradation rate (k MD ) as an independent model parameter (equal to, greater or less than the endogenous degradation rate, respectively, in the three scenarios described above). This basic mechanistic model can be further customized to include, e.g., competition with metabolites or endogenous ligands (Fig. 3) -which may be relevant for early PROTACs whose PoI ligand consists of a pre-existing small molecule inhibitor which binds to an active site of the target protein [22], or by unfolding the PROTAC degradation rate in a series of transit compartments to better describe the ubiquitination or de-ubiquitination process, as done, e.g., in [23,24]. Since the relevance of these model components may be target or chemo-type specific, for the sake of generality they will not be included in this analysis, although the methodology utilized here can be applied to an extended version of the model as well (provided the assumption of total PROTAC and total ligase conservation is met).
The governing equations are derived from mass balancing principles and they read as follows: Fig. 2 Ternary complex levels as a function of PROTAC concentration in a pure binding system (left); target levels relative to baseline as a function of PROTAC concentration in a biological setting (right) with initial conditions Tð0Þ ¼ T 0 , Lð0Þ ¼ L 0 , Pð0Þ ¼ P 0 , PLð0Þ ¼ PTð0Þ ¼ TPLð0Þ ¼ 0, where T, L, P, PT, PL, TPL stand for target, ligase, PROTAC, PROTAC-target complex, PROTAC-ligase complex, and ternary complex concentration, respectively, and will be referred to as states of the system (dependence on time t has been suppressed to ease the notation). Model parameters are defined in Table 1. Note that on/off rates are correlated via cooperativity a as [25]: and since the proportionality constant a is the same in (5), the relationship between on/off rates can be synthetically expressed as This means that, even though the binding kinetics is governed by 8 parameters, only 7 of them are independent. In other words, knowing any 7 on/off rates is sufficient to calculate the remaining one via Eq (6). Adding endogenous synthesis and degradation rates (k syn , k deg ) as well as binary   and ternary complex degradation rates (k MD , k PRO ) gives a total of 11 independent parameters. Under the assumption of total PROTAC and total E3 ligase conservation the differential equations for free PROTAC and free E3 ligase are redundant as they can be obtained from conservation laws as where L 0 and P 0 are the ligase baseline and PROTAC concentration, respectively. The steady state of the system is obtained by setting each derivative in (4) to 0, and the same method described in [18] previously adopted for MDs can be applied here.

Global sensitivity analysis
Sensitivity analysis is a powerful tool to understand how and to what extent each model parameter (input) affects the response (output) [13]. Local sensitivity analysis studies the system states variability as a single parameter varies, all the others being fixed. While this approach is extremely convenient for its simplicity of implementation and easiness of interpretation, it can be misleading as the sensitivity of the response to one parameter can depend on the values of all the other fixed parameters. In other words, the model output can be sensitive to a parameter p H for a given set of the other parameter values and at the same time insensitive to p H for a different set of fixed values. Therefore, this approach can be useful and unbiased only if confidence in the fixed parameters is high.
Differently, global sensitivity analysis studies the output variability over the whole parameter space, i.e. as all model parameters are changed simultaneously. As a result, the response variability characterization is more robust as it only depends on the assumed or observed distribution of each parameter on a given feasible range rather than on a single value [26]. At the same time, though, visualizing the response variability and quantifying the contribution of each parameter to it is increasingly challenging with the dimension of the parameter space: it is well known that the number of points to accurately sample a parameter space grows exponentially with the dimension of the parameter space itself (''curse of dimensionality''). In other words, if N samples are sufficient to describe the distribution of a single parameter, an order of N P points will be required to equivalently accurately sample P parameters. For instance, if 10 samples are used for each parameter of the PROTAC model (4) the total number of required samples would be approximately 100 billions (10 11 ).
Each set of the 100 billion parameter combinations will generate a model output, and visualizing and interpreting 100 billion model simulations is clearly more challenging than plotting 10 of them (as a single parameter changes), as well as more computationally expensive.
A plethora of mathematical tools to quantify the impact of each parameter on the response and to tackle the curse of dimensionality are available. In this work Sobol indices [27] are used to assess the fraction of total variability associated with each parameter, which is a random variable represented via Polynomial Chaos Expansions (PCEs) [28][29][30][31]. Parameters are assumed to be uniformly distributed around a given mean and standard deviation, and are hence accordingly represented by first-order PCE of Legendre polynomials, which maximize the convergence rate according to the Askey scheme [32,33]. Parameter uncertainty is propagated in the system via Non-Intrusive Spectral Projection (NISP) [34]. As a result, the stochastic model output can be represented as a PCE as well, whose coefficients can be easily used to calculate Sobol indices. In order to reduce the computational cost of sampling associated with NISP without compromising accuracy, Smolyak sparse quadratures are employed [35]. The sparsity level has been manually increased until no significant change in the Sobol indices estimates was observed. Uncertainty propagation via NISP and Sobol indices calculation was handled with the C?? library UQTk developed at the Sandia National Laboratories [36], embedded in a MATLAB implementation of the mechanistic PROTAC model (4).

In vitro
Test compounds were evaluated at 12 concentrations obtained with a 1:3 dilution factor, with two replicates per concentration. Remaining PoI levels were measured via Western Blots, ERD9 [37], Immunofluorescence [38] or HiBit technology [39] and expressed as a fraction of protein in DMSO treated cells (i.e. baseline).
The following bi-sigmoidal model describing the remaining fraction of target at steady state ( b T SS ) as a function of concentration (C) was fitted to each individual dose-response at steady state to capture any potential hook effect and calculate maximal degradation and potency: b Eq (8) describes the response as a combination of two sigmoidal curves with half-maximal concentrations EC 50 and IC 50 , respectively. E max represents the overall maximal effect, while E loss can be interpreted as the fraction of degradation lost to the hook effect. To avoid overparametrization the Hill coefficient of the sigmoid corresponding to the hook effect was fixed to 1, while it was estimated (h) for the sigmoid corresponding to increasing degradation. Because the response is the result of the contribution of two distinct sigmoidal curves, the observed potency (DC 50 ) may not exactly correspond to EC 50 , therefore it was calculated numerically from (8) as the concentration delivering half-maximal effect. The concentration corresponding to maximal degradation (DC max ) was also calculated numerically as the root of the first derivative, and maximal degradation as Whenever dose-response time courses were available, model (8) was embedded into the following turnover model describing the fraction of target at any given time ( b T): where the endogenous fractional turnover rate k deg was estimated or fixed to experimental data obtained from Stable Isotope Labeling with Amino acids in Cell culture (SILAC) [40]. Experimental data in Sect. 4.1.3 has been generated with AstraZeneca proprietary Selective Estrogen Receptor Degraders (SERDs) [41][42][43][44][45].
Endogenous PoI or E3 ligase levels in different cell lines were assessed via Western Blots.

In vivo
In vivo data was generated in NSG mice implanted with a patient-derived xenograft tumor model. C-PROTAC-006 was dosed orally on a daily schedule for three days at 30 mg/ kg, 60 mg/kg or 100 mg/kg. Plasma concentration at different time points was quantified via Liquid Chromatography with tandem Mass Spectrometry (LC-MS/MS). Remaining protein levels relative to vehicle baseline levels were assessed by Western Blots at 6h, 24h, 48h after the last dose.
A one-compartmental pharmacokinetic model with firstorder absorption was fitted to the plasma concentration data and used as a driver of the pharmacodynamic model (9) parametrized from in vitro data generated in the same cell line to obtain predictions of in vivo degradation kinetics.

Exact steady state solution
By applying the method in [18] the following expression can be obtained, which implicitly describes the free target at steady state (see Appendix B.1): Since total protein amounts (T þ TD) relative to baseline are more easily accessible experimentally, through some manipulations an expression for the fraction of total remaining target can also be obtained as a function of the free MD concentration (see Appendix B.2): Notice how Eq. (11) has the structure of a decreasing sigmoidal curve of the form with top plateau and Hill coefficient equal to 1, and bottom plateau ( b T min ), maximal degradation ( b D max ) and concentration for half-maximal degradation (DC 50 ) respectively equal to which is consistent with [46].

Exact steady state solution
The steady state solution for PROTACs is convoluted due to the complexity associated with three-body binding kinetics. In fact free target T and free ligase L are obtained simultaneously by numerically solving the following non-linear system of conservation laws for PROTAC and E3 ligase 1 : where P, PL, PT and TPL are all known calculated functions of T and L encapsulating surrogates of the original mechanistic parameters (Appendix D.1). Solving (14) means finding the one pair T H ; L H À Á among all possible system states that does not violate the conservation laws, i.e., if we call Z 1 ðT; LÞ ¼ L tot ðT; LÞ À L 0 and Z 2 ðT; LÞ ¼ P tot ðT; LÞ À P 0 the amount by which each conservation law is violated, the solution to (14) will satisfy Z 1 ðT H ; L H Þ ¼ 0 and Z 2 ðT H ; L H Þ ¼ 0 simultaneously. Figure 4 provides a visual example: Z 1 and Z 2 are surfaces whose shape is dictated by the mechanistic parameters; the intersection between the two surfaces and the horizontal plane Z ¼ 0 corresponds to the pair T H ; L H À Á that satisfies (14). Once T H and L H are calculated numerically, all the remaining states as well as the total remaining target which is typically available through measurements, can be easily obtained through the exact expressions (D25a), (D25b), (D25c) (Appendix D.1). In particular, the total remaining target as a fraction of the baseline reads as: and it can be verified that it is equivalent to the numerical steady state solution of the differential Eq. (4) (Fig. 4).

Global sensitivity analysis
Global sensitivity analysis was performed on the mechanistic model (4) assuming uniform distribution of the parameters over the ranges listed in Table 2. In order to reduce the exponentially high computational cost of GSA on a high-dimensional parameter space, no cooperativity was assumed (a ¼ 1, or k PL Smolyak sparse quadratures of level 6 were used, corresponding to 242815 sampled parameter combinations. Sobol indices of DC 50 and D max were calculated for each parameter from the corresponding simulated doseresponses and grouped by category to ease the interpretation (Fig. 5). More precisely, the fractions of variability due to binding kinetics, degradation kinetics, and target baseline sum up the contribution of the Sobol indices of on/ off rates (k T on , k T off , k L on , k L off ), degradation rates (k MD , k PRO ), and endogenous target synthesis/turnover (k syn , k deg ), respectively. Most of the variability of both DC 50 and D max is due to changes in ligase baseline ( [ 50%), followed by degradation efficiency ( $ 20 À 40%) and target baseline ( $ 20%). Based on the selected parameter range and distribution, binding kinetics contributes only minimally to the overall response variability ( $ 5%).

Compound optimization
Binding equilibrium constants (K d ¼ k T off =k T on ) are often used to drive SAR of small molecules. While they have proved useful in practice [47][48][49][50], from a mechanistic viewpoint they can only provide limited understanding of the major drivers of pharmacology: 1. On/off rates matter individually. Consider the case where k T off is minimized vs. the case where k T on is maximized. In both cases the overall effect on thermodynamics is that the binding affinity increases (i.e. K d decreases), however the same variation in affinity can actually result in different degradation profiles. This is particularly evident in the two extreme scenarios where k T off is extremely low (namely k T off ! 0 þ , as it could be the case, e.g., of a covalent binder) or when k T on is much greater than k T off þ k MD (namely k T on ! þ1). The definition in (13) clearly shows that by reducing k T off one can only lower DC 50 to ð1 À b D max Þk MD =k T on , whereas by increasing k T on there is virtually no limit to how potent a compound can be made (namely, up to DC 50 ¼ 0). Analogously for b D max , from Eq. (11) we can clearly see that, for the same degradation efficiency (k MD ), ever faster on-rates will reduce the target to a minimum relative amount equivalent to the ratio of endogenous and MD degradation rate (Eq. (17a)); differently, ever slower off-rates will spare a portion of residual protein that is inversely proportional to the amount of free drug at steady state (Eq. (17b)): In formulae, notice that the limit T OFF in (17b) is equal to T ON (17a) added by a residual term, which clearly shows that optimizing off rates leads to inferior degradation (Fig. 6). 2. Maximizing on-rates is the most efficient way to maximize binary complex degradation via binding kinetics. Minimizing (11) solely via binding kinetics entails minimizing the residual term. Because of its functional structure, this can only be achieved by making the ratio k T off þ k MD k T on as close to 0 as possible, i.e., by maximizing k T on . 3. Binding and degradation optimization is coupled.
Eq. (11) highlights how the residual term that prevents a MD from achieving maximal degradation (1 À k deg =k MD ) is governed not only by the binding equilibrium constant k T off =k T on but also by the ratio of binary complex degradation rate to the binding on-rate (k MD =k T on ), which parameters are clearly different in nature (pharmacology vs. chemistry). This shows that binding and degradation cannot be thought of separately. In other words, relatively slow binding can be efficient enough for a slow degrader, but on the other hand a fast degrader will require fast binding to deliver its full potential. For many years medicinal chemist have focused on increasing the drug-target residence time by means of optimizing off-rates [47,48]. While that approach has been successful, this analysis suggests that for degraders optimizing on-rates is a more efficient approach. Copeland et al already highlighted that optimization of the target residence time (i.e. off-rate) may have limited utility in cases where the rate of new protein synthesis (and degradation) plays an important role [47,48]. Multiple authors [47,48,51] mention an upper limit to on-rates (10 8 À 10 9 M À1 S À1 ) given by the rate of diffusion of the two binding partners in physiological solutions, however Schoop et al [51] show that several discovery programs have not yet achieved such upper limit, indicating that there may still be room to further increase on-rates. Medicinal chemists have been confronted with the question of whether it is possible to design on-rates for a given target within a given compound series but Schoop has shown different compound series of distinct targets which display significant differences in their on-rates. Obviously the difficulty facing medicinal chemist is establishing general k on optimization strategies. Optimization of on-rates through SAR has been limited and mostly empirical due to the difficulty to predict the physicochemical steps involved in receptor-ligand association (protein conformational re-arrangements, protein desolvation and molecular orbital orientation) [47,48,51]. 4. Improving degradation efficiency is necessary to increase maximal degradation, and sufficient but not necessary to improve degradation potency Because D max only depends on the (endogenous and) MD degradation rate via (13), the only way to increase maximal degradation is by increasing k MD . Moreover, because degradation potency directly depends on D max (13), improving degradation efficiency (k MD ) always results in higher potency (lower DC 50 ) as well. However, the opposite is not necessarily true: an improvement in potency driven by optimized binding kinetics (faster k T on or slower k T off ) has no impact whatsoever on D max . 5. Complete degradation can only be reached asymptotically; faster turnover proteins are harder targets. It is impossible to achieve net 100% degradation because b T min ¼ k deg =k MD is always strictly positive, no matter how small the ratio. Moreover, faster protein turnover requires higher degradation efficiency to maintain a given level of target knock-down.
In conclusion, optimizing binding affinities is a simplistic approach that can lead to ambiguous outcomes because the response is sensitive to on and off rates individually, and binding kinetics should be optimized relatively to degradation kinetics (and conversely).

Model identifiability and data requirements
Although analyzing the fraction of total remaining target relative to baseline would be more relevant because truly representative of the available experimental data, for the sake of simplicity we now focus on the expression of the free target concentration. In Appendix B.3 we show that the same conclusions hold true for the fraction of total remaining target.
If we lump the coefficients in (10) into a new synthetic parametrization to highlight the equation structure (Eqs. (18)-(19)) it becomes evident that while the original mechanistic model (1) is governed by five parameters (k syn , k deg , k T on , k T off , k MD ), the corresponding steady state (10) only features three independent degrees of freedom (a, b, c, see Appendix B.3), which are surrogates of the original mechanistic parameters: This apparently simple observation is extremely informative on model identifiability and data requirements: 1. If the only available data is free target concentration at steady state T, the surrogate parameters a, b and c can be uniquely estimated, but none of the original mechanistic parameters can be uniquely identified. 2. If, in addition to point (1), the endogenous degradation rate k deg is available, the MD degradation rate k MD can be obtained from (19) as and consequently the endogenous synthesis rate k syn and protein baseline T 0 as 3. If, in addition to points (1) and (2), the binding equilibrium constant K d ¼ k T off =k T on is available, also the individual on/off rates can be calculated as Example: baseline estimation from degradation data Eq. (B14) in Appendix B.3 describing the fraction of total remaining target at steady state was fitted in Matlab to the in vitro degradation dose-response of 12 SERDs to obtain the surrogate parameters a, b, c (An example is shown in Fig. 7, left). The measured endogenous fractional turnover rate of ERa in MCF-7 (k deg ¼ 0:15/h, Appendix C) was used to calculate the MD degradation rate k MD via (20). Then, ER baseline was estimated from each compound from (21) with an interquartile range of 1:86 À 4:51nM (Fig. 7, right), which is within published concentrations of ER for MCF-7 cells [52,53].

Compound optimization
We have shown that optimizing MDs based on binding affinity alone is sub-optimal because the same change can produce different effects on degradation, depending on whether such variation is due to slower/faster on-rate vs off-rate. Differences in the effect can be traced back to binding kinetics and degradation kinetics being coupled or, more precisely, to the proportion of MD degradation rate k MD and on-rate k T on . In this respect, although drawing similar conclusions for PROTACs is less straightforward due to the complexity of the steady state solution, understanding parameter dependencies (i.e., proportions of mechanistic parameters or which ones matter with respect to which others) can provide at least general guidelines on the most critical aspects of optimization.
The core of PROTACs inter-parameter relationships and proportions is encoded in coefficients c ij 's and t ij 's (Appendix D.2). Rather than analyzing the expression of each individual coefficient, we will summarize such relationships in a table where both rows and columns correspond to kinetics, pharmacological or biological parameters, and the cell in a given row r, column c is filled if the proportion between the two parameters in the corresponding row and column has an effect on steady state degradation (i.e., if they appear in a ratio in the coefficients). For instance, the MD solution (10) contains the ratios k T off =k T on , k deg =k MD , k MD =k T on , which correspond to the table pattern in Fig. 8 (top). The non-empty fingerprint of the top-right or bottom-left offdiagonal blocks clearly marks the coupling between binding and degradation kinetics previously noted (k MD =k T on ). Furthermore, while the top-left diagonal block encodes the binding affinity k T off =k T on , the bottom-right one captures the maximum achievable degradation (17a), given by the proportion between the endogenous and MD degradation rates k deg =k MD , thus coupling biology with pharmacology.
The fingerprint of the PROTACs solution is shown in Fig. 8 (bottom). As for MDs, the table is partitioned into binding and degradation kinetics, however the former block is further sub-divided into target and E3 ligase kinetics of single species or binary complexes, for a total of 2 Â 2 sub-blocks of size 4 Â 4.

Binding and degradation optimization is coupled.
First and foremost, as in the previous example it is evident that the matrix pattern is not block-diagonal, i.e., binding and degradation kinetics are intertwined. 2. PROTACs are more than just the sum of their components. More specifically, even the top-left two-bytwo target/ligase binding kinetics block is not blockdiagonal: this means that PROTACs-induced response is the result of not only target and ligase kinetics individually, but also of the interplay between the two. For instance, target-PROTAC on-rate (k T on ) matters relatively to ligase-PROTAC on rate (k L on ). This fact underscores that PROTACs design goes way beyond the optimization of its individual components (E3 ligase and target warheads), i.e., optimizing each portion independently may not guarantee an optimal compound overall.
3. On/off rates matter individually. Each on/off rate plays a unique and specific role in relationship with the rest of the parameters (the fingerprint of each row or column is specific to each parameter). As for MDs, even more so for PROTACs this implies that changes not only to each binding affinity, but also on cooperativity can produce a different effect depending on whether on-vs off-rates are changed. Therefore studying how the response varies as cooperativity changes is an ill-posed question as it is for binding affinities. 4. Increased cooperativity via faster on-rates promotes degradation. As described above, an excess of PT (or PL) binary complex can impair efficient degradation up to causing protein stabilization. The exact solution fingerprint of k MD shows that PT accumulation can be mitigated via binding kinetics in two ways: by favoring the ligase pathway via faster PROTAC-ligase on rates (k MD =k L on ratio), or by increasing cooperativity via on rates (k MD =k PL on , k MD =k PT on ratios).

Model observability
Despite the complexity of calculations and resulting expressions, the compact notation in Eq. (14) is key to highlight a structural feature of the system with implications on minimal data requirements. First of all, the E3 ligase baseline L 0 is an input to the system, just as much as the drug concentration P 0 , therefore experimental knowledge of it is a prerequisite. Secondly, if we first consider the simpler case of MDs, expression (18) clearly shows that observations of only one species (the remaining free target T following a corresponding given dose D 0 ) is necessary and sufficient to infer the remaining states of the system (binary complex TD and free MD D from (B4a) and (B5) in Appendix B, respectively) and possibly estimate the surrogate parameters (provided enough distinct observations). Differently, for PROTACs, Eqs. (14) feature two independent variables (T and L) for each given concentration P 0 which cannot be further reduced. This suggests that observations of at least two different states are required in order to infer the state of the whole system and possibly estimate its surrogate parameters. Such two species do not necessarily need to be free target and free ligase: any other pair of observations (except for total ligase and total PROTAC, which are conserved) is viable because T and L can be retrieved from relationships (D23a), (D25a)-(D25c) (Appendix D.1) to then solve (14). For instance, if the concentrations T H and T H tot were observed for free target T and total target T tot , free ligase concentration L could be calculated by solving T H þ PTðL; T H Þþ TPLðL; T H Þ ¼ T H tot , where the notation ðL; T H Þ indicates functional dependence on ligase, given the target concentration, and the exact functional relationships of PT and TPL are known from (D23a) and (D25c). Analogously, if total target T H tot and ternary complex TPL H were observed, target and ligase concentrations could be retrieved by solving Eq. (15) and (D23a) for T and L simultaneously, given T H tot and TPL H .
In summary, the exact solution structure practically implies that measuring total remaining target concentration only is insufficient to infer the state of the whole steady state system, which also requires knowledge of the E3 ligase baseline.
Model identifiability While these conclusions characterize the steady state system observability, i.e., the ability of inferring the state of the whole system given a set of observations, understanding which surrogate parameters can be uniquely identified (such as a, b and c for MDs) is not as immediate. Ultimately, the steady state solution is characterized by 14 Fig. 8 Tables of parameters interactions emerging from the exact steady state solution of MD (top) and PROTAC (bottom) mechanistic models. A filled cell in row r, column c indicates that the parameters in row r and column c matter relatively to each other. Relationships marked with an asterisk refer to an equivalent dual parametrization of the exact solution (not shown, see Appendix D) and are included for completeness surrogate coefficients (c ij 's and t ij 's, Appendix D.2) which, alone, are insufficient to retrieve the original 11 mechanistic parameters. In fact, by constructing a Gröbner basis [54] it can be shown that only 9 of them are actually independent. This means that at least 2 mechanistic parameters (or corresponding surrogates) have to be measured to be able to retrieve all the others from 9 model coefficient estimates. For instance, if the endogenous degradation rate k deg (or the endogenous synthesis rate k syn , or the binary complex degradation rate k MD ) and the PROTAC-target binding affinity k T on =k T off are available from measurements, all the original mechanistic parameters can be calculated (Appendix D.3). The generalization of this conclusion to all possible sets of measurements from which it is possible to retrieve the mechanistic parametrization will be addressed in future work.

Data prioritization: impact of protein endogenous baselines on pharmacology
While drawing firm guidelines on PROTACs optimization is challenging due to the complexity of the system, the outcome of global sensitivity analysis can be extremely valuable to guide data prioritization and model simplification. First and foremost, it is striking how most of the total response variability (up to 75%) is not driven by chemical or pharmacological factors such as binding and degradation rates, which can indeed be optimized, rather by biological factors, i.e. ligase and target baselines, over which we have no control. This shows how critical it is to characterize the distribution of such parameters in the patient population and, consequently, to judiciously select the most appropriate pre-clinical patient representative tumour models [55]. More precisely, data collected from four different AstraZeneca PROTACs projects clearly show how D max correlates with the ratio of ligase and target baselines, i.e., the higher the ligase levels relatively to the target the higher the degradation (Fig. 9). It is to be noted that here baseline quantification is obtained from Western Blots and hence is only semi-quantitative. While efforts are in place to generate absolute quantification of protein levels in relevant models to confirm such relationship, the fact that a similar trend is observed across multiple targets and compounds from different chemical series is early evidence that such correlation may hold true in practice, both in vitro and in vivo [15]. Note that even within the same target space not all molecules may display the same type of correlation: a range of different slopes or D max may be observed. For instance, Fig. 9 (center) shows how A-PROTAC-001 displays a flatter profile with a lower D max compared to other compounds targeting the same protein of interest and ligase. What makes a PROTAC more or less sensitive to changes in baseline ratio is still unclear, and which profile is more desirable depends on the PKPD/E and toxicity relationship of the target. In this example, if 80 À 90% maximal degradation is sufficient to drive efficacy then A-PROTAC-001 is preferable because it is expected to be efficacious in most of the relevant cell lines and hence ultimately in a wider set of baseline ratios (which may be observed in the target population). On the other hand, if nearly complete degradation is required for efficacy then A-PROTAC-002-005 may be preferable because they can achieve higher D max up to $ 100%. However, in this scenario characterizing the ligase and target baseline distribution in the patient population is even more critical: if the ratio is relatively low in most patients these PROTACs may be unable to achieve complete degradation due to the steep correlation with the baseline Fig. 9 Correlation between D max and ratio of ligase and target baselines across different targets (left; only one compound per program is shown for easiness). D max /baseline ratio correlation of several compounds targeting protein A (center). DC 50 /baseline ratio correlation of several compounds targeting protein C (left) ratio, and as a result efficacy may be compromised. At the same time, though, a sharp slope may enable a better therapeutic index: if the relevant tox tissues display lower ratios compared to the patient representative tumour models a steeper relationship would entail degrading the target in the tumor while sparing it in normal tissues [15].

A pragmatic approach to modelling
This analysis implies an empirical approach to build a mathematical model of medium complexity with the potential of being truly predictive. While, on the one hand, the fully mechanistic model is over-parametrized and hence unlikely to be predictive, on the other hand traditional turnover models such as (9) are also unlikely to be predictive across patient representative tumour models because unable to capture differences in D max driven by changes in baselines. However, GSA of the former suggests that it is possible to ''enrich'' or ''augment'' the latter with data-driven components. More precisely, D max can be implemented as a (sigmoidal) function of the baseline ratio as inferred from the data rather than as a constant parameter, and the same can be done for DC 50 (Fig. 9, right). For instance, if q is the ligase/target baseline ratio these two parameters can be described as DC 50 ðqÞ ¼ðDC 50;max À DC 50;min Þe ÀjÁq þ DC 50;min ð23Þ and used in Eq. (9). Clearly, this type of modelling requires generating dose-responses in multiple patient representative tumour models spanning a range of baseline ratios, which is more demanding than a simple dose-response but still less prohibitive and more easily accessible than individual on/off rates. Not only, it has the potential of predicting the response of new patient representative tumour models or tissues (or even patients) as they become of interest or available.
Whenever baseline and/or dose-response data is unavailable or insufficient to fully parametrize relationships (23), it is critical to use the same in vitro cell line as the in vivo model to increase the chances of operating in biological settings with comparable baseline ratio q (even without knowing its value) while mitigating the risk of a disconnect in D max and potency due to potential differences in baseline ratio. For instance, Fig. 10 shows how a PD model built on degradation time course data in model C-9 in vitro, driven by the plasma PK profile fitted in vivo in the same C-9 model, can successfully predict degradation kinetics in vivo at multiple dose levels.

Conclusion
In this work we have shown the value of an integrated modelling approach for degraders which combines the benefits of traditional turnover models and fully mechanistic models to address two key project questions in drug discovery programs, i.e. (i) how to drive compound optimization and (ii) how to predict pharmacology across patient representative tumour models or the kinetics of in vivo degradation from in vitro data, ultimately to support dose selection and predictions.
Whenever an exact mechanistic steady state solution is available it can be utilized to precisely understand the role of each parameter on the response and hence which kinetic ''knobs'' (e.g. on/off rates) need to be tuned to achieve a desired pharmacological profile (D max , DC 50 ). If such properties can be measured experimentally in highthroughput format they can directly inform and guide the design of novel compounds, otherwise it may be possible to obtain them indirectly from degradation data.
Following this approach we have shown how on/off and degradation rates are related to potency and maximal effect of monovalent degraders, and how such relationship can be used to suggest a compound optimization strategy based on faster association rather than slower dissociation. Moreover, the same relationship indicates which additional experimental data (e.g. endogenous protein turnover and binding affinity) should be collected to retrieve all the original mechanistic parameters, which we demonstrated for SERDs. Whenever an exact steady state solution is not available, or when its complexity is prohibitive (such as for PROTACs), an empirical data-driven, mechanism-agnostic approach offers a sub-optimal yet more pragmatic option. Nevertheless, even complex exact steady state solutions can provide insight on the type of observations required to ensure the predictive capacity of a mechanistic approach. Specifically for PRO-TACs, the structure of the exact steady state solution suggests that the total remaining target at steady state, which is easily accessible experimentally, is insufficient to reconstruct the state of the whole system at equilibrium and observations on different states (such as binary/ternary complexes) are necessary. Not only, it highlights how the E3 ligase baseline levels are a direct input to the system -just as the compound concentration -and therefore should be measured beforehand in biological models of interest.
While parameter identifiability is challenging for PRO-TACs, global sensitivity analysis suggests that both target and ligase baselines (actually, their ratio) are the major sources of variability in the response of non-cooperative systems, which speaks to the importance of not only generating such data as early as possible in drug discovery, but also characterizing their distribution in the target patient population to maximize the probability of response in the clinic. Importantly, GSA can also suggest pathways to enrich (too) simplistic turnover models just enough to unlock their predictive capacity. The pragmatic approach we proposed here is being applied within AstraZeneca to PROTAC programs, accelerating the drug discovery pipeline and increasing the chances of success in the clinic.
The first applicability requirement for method [18] is N P N C þ N B , which is met (2 1 þ 1). Although the total number of Equations is N E ¼ 3, the linear conservation law introduces constraints among states that reduce the number of independent equations to 2 (as shown in Eq. (3)). Therefore we will solve for N D ¼ 2 dependent variables and N F ¼ N V À N D ¼ 2 free variables.
Secondly, we need to find a clean set of DN ¼ N P À N C ¼ 1 variable that appears in bilinears in such a way that each bilinear term contains at most one of them. Because there is only one bilinear term, T Á D, T or D are equivalent options. For example, set X X ¼ fTg. The clean set is key to the applicability of the method, as it enables the linearization of the non-linear component of the system.
It is now possible to identify a clean partition of free and dependent augmented variables. Dependent variables (X D ) will include states that do not appear in bilinears (X M ) and the bilinear term. Free variables (X F ) will include states that appear in bilinears (X P ) (in this case there are no remaining bilinears since the only one has been included as a dependent variable). In summary: With such partition in place, (B1) can be rearranged as a system of 2 independent linear equations in 2 unknown dependent variables X D parametrized on the free variables X F (on the right-hand side): Solving such linear system provides expressions for the binary complex TD and the bilinear term as functions of the remaining states T and D: Because the bilinear term contains only one variable of the clean set X X (i.e., it is linear in the clean set), D can be easily obtained as function of T: At this point, by plugging in expressions (B4a) and (B5) the conservation law (2) can be expressed by and solved for for T only: which is the implicit exact solution (10). In summary, this method linearizes the original nonlinear system (B1) by treating bilinear terms as additional (dummy) variables and by wisely choosing a subset of them as parameters so that the non-linear core of the system is curtailed to the sole conservation laws.

B.2 Exact expressions of total target relative to baseline as a function of free drug
Most in vitro assays detect the (fraction of) total amount of remaining target T tot ¼ T þ TD, whose exact expression as a function of T can be easily obtained from (B4a): Kinetics parameters such as k T on and k T off do not appear explicitly, rather implicitly through T. Although finding a formula for T is impossible -it can only be calculated numerically from (B6) -relationship (B4b) can be used to express T in terms of D, whereby on and off rates appear explicitly: Plugging (B8) into (B7) and dividing by the baseline T 0 ¼ k syn =k deg we obtain the fraction of total remaining target relative to baseline: which is Eq. (11).

B.3 Minimal parametrization of total remaining target relative to baseline
Since the available in vitro data is typically in the form of total target as a fraction of baseline, in order to enable parameter estimation we need to express the exact solution (B6) in terms of b T tot ¼ T tot =T 0 . This task can be performed in three steps: 1. Express (B6) in terms of the normalized free target b T ¼ T=T 0 . By multiplying and dividing each term containing T in (B6) by the baseline T 0 ¼ k syn =k deg (note that this equates to multiplying by 1, which does not change the equation) we obtain: and recalling the definitions of the surrogate parameters a, b, c (19): 3. Plug (B13) into (B11). We obtain which can be expanded as a polynomial function in b T tot : The surrogate parameters a, b, c are still identifiable from (B15), in fact they can be uniquely obtained by estimating which are certainly identifiable by construction, and by setting Appendix C SILAC modelling The following mono-exponential model was fitted to heavy-and light-labelled SILAC data in untreated MCF-7 cells (Fig. 11), with the endogenous fractional turnover rate k deg as a shared parameter: Model estimates are reported in Table 4 with 95% confidence intervals.

D.1 Steady state solution
The steady state of the dynamical system (4) is obtained by setting the time derivatives to 0: Let us define equations and variable sets as in Table 5, with corresponding cardinality. The first applicability requirement on dimensionality N P N C þ N B is met (5 2 þ 4). Although the total number of equations is N E ¼ 6, the two linear conservation laws introduce constraints among states that reduce the number of independent equations to 4. Therefore we will solve for N D ¼ 4 dependent variables (X D ) and N F ¼ N V À N D ¼ 6 free variables (X F ). In order to be able to apply the method, the choice of free and dependent states cannot be fully arbitrary. We need to find a clean set of DN ¼ N P À N C ¼ 3 variables that appear in bilinears in such a way that each bilinear term contains at most one of them. For instance, P, PT and PL all appear in the bilinears X B , but each bilinear contains each one of them only once (in other words, P, PT and PL are never multiplied by one another). Therefore, X X ¼ fP; PT; PLg is a clean set for system (D19). The clean set is key to the applicability of the method, as it enables the linearization of the non-linear component of the system.
It is now possible to identify a clean partition of free and dependent augmented variables. Dependent variables will include states that do not appear in bilinears (X M ) and a subset of DN ¼ 3 bilinears such that each element of the clean set X X appears once and only once in all the selected bilinears, for a total of N D ¼ 4. For instance, there is one and only one bilinear in the set X B;X X ¼ fL Á P; T Á PL; L Á PTg that contains each one of P, PT and PL (i.e., elements of the clean set X X ) once: only L Á P contains P, only T Á PL contains PL, and only L Á PT contains PT. 2 Free variables will include states that appear in bilinears (X P ) and the remaining bilinears (X B n X B;X X ). In summary: With such partition in place, system (D19) can be rearranged as a system of 4 independent linear equations in 4 unknown dependent variables X D parametrized on the free variables X F (on the right-hand side). Solving such linear system provides expressions for TPL (i.e., X M ) and the bilinears X B;X X as functions of the remaining states T, P, L, PT, PL (i.e., X P ) and the bilinear T Á P (i.e., X B n X B;X X ). Thanks to the rational choice of free and dependent variables, the three equations associated with the bilinears X B;X X in turn make up a linear sub-system in the clean variables of X X , parametrized on T and L (i.e., X P nX X ). As a result, solving this linear sub-system provides expressions for P, PT and PL in terms of T and L, which in turn allows for reformulating the conservation laws as functions solely of the two latter variables. Finally, solving the conservation laws for T and L allows the retrieval of all the other states (and dummies). In summary, this method linearizes the original nonlinear system by treating bilinear terms as additional (dummy) variables and by wisely choosing a subset of them as parameters so that the non-linear core of the system is curtailed to the sole conservation laws.  Isolating dependent variables on the left-hand side and free variables on the right-hand side in (D19) yields: where the equations for P and L have been dropped due to redundancy (they can be retrieved from the remaining states and the conservation laws) and each bilinear term has been highlighted in brackets as a whole dummy variable. System (D21) can be written in matrix form as Eqs. (D23b)-(D23d) form a linear system in the clean set X X ¼ fP; PT; PLg, parametrized on T and L, which can be recast in matrix form as from which T and L can be numerically calculated for any given concentration P 0 and ligase baseline L 0 .

D.2 Steady state solution coefficients
The rationale for coefficients notation is displayed in Table 6. The subscripts in each coefficients are determined as follows: • c ij is the coefficient in row i, column j of the clean set system in matrix form (D24); • t ij , p ij , b T;ij , b L;ij , d ij multiply the monomial term T i Á L j in the corresponding formula as in Table 6; • t H multiplies the PT binary complex in (D23a) (and is the only coefficient that does not multiply a monomial term in two variables of the form T i Á L j ).

D.3 Mechanistic parameters retrieval
Assume endogenous degradation rate k deg and the PRO-TAC-target binding affinity k T on =k T off are known form measurements, and that the 9 coefficients t 00 , t 10 , t H , c 11 , c 13 , c 20 , c 23 , c 30 , c 33 are estimated. Then, synthesis and degradation rates are obtained as: From the definitions of c 11 and c 30 , and given the ratio k T on =k T off , it follows k T on and finally Acknowledgements The authors would like to acknowledge multiple AstraZeneca Targeted Protein Degradation discovery teams for the generation of both in vitro and in vivo data in this work.
Author contributions All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by S.G. The first draft of the manuscript was written by S.G. and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Declarations
Conflicts of interest SG and PMG are AstraZeneca employees and hold shares and stock options.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.