A Kinetic Model Explains Why Shorter and Less Affine Enzyme-recruiting Oligonucleotides Can Be More Potent.

Antisense oligonucleotides complementary to RNA targets promise generality and ease of drug design. The first systemically administered antisense drug was recently approved for treatment and others are in clinical development. Chemical modifications that increase the hybridization affinity of oligonucleotides are reasoned to confer higher potency, i.e., modified oligonucleotides can be dosed at lower concentrations to achieve the same effect. Surprisingly, shorter and less affine oligonucleotides sometimes display increased potency. To explain this apparent contradiction, increased uptake or decreased propensity to form structures have been suggested as possible mechanisms. Here, we provide an alternative explanation that invokes only the kinetics behind oligonucleotide-mediated cleavage of RNA targets. A model based on the law of mass action predicts, and experiments support, the existence of an optimal binding affinity. Exaggerated affinity, and not length per se, is detrimental to potency. This finding clarifies how to optimally apply high-affinity modifications in the discovery of potent antisense oligonucleotide drugs.

Introduction

Antisense oligonucleotides can be used to elicit RNase H-mediated cleavage and decay of target RNA.[1] RNase H is a non-sequence-specific enzyme that recognizes RNA–DNA heteroduplexes and specifically cleaves the RNA strand.[2] Although the basic principle behind antisense was realized early,[3] the binding affinity and pharmacokinetics of natural, unmodified, oligonucleotides were found to be insufficient for their use as systemic drugs.[4] Therefore, central to the development of oligonucleotides for therapeutics lies efforts in medicinal chemistry to improve oligonucleotide stability, biodistribution, as well as RNA binding affinity. With the introduction of high-affinity locked nucleic acids (LNAs),[5,6] oligonucleotides as short as 12nt in length have been reported to achieve sufficient binding affinity to potently silence their targets.[7,8] Other types of modifications, such as 2′-O-methoxyethyl and 2′-O-methyl,[9] are also in use, although here 20nt are usually needed for sufficient binding affinity. Since none of these modifications allow recruitment of RNAse H, modified oligonucleotides are typically designed as gapmers, with high-affinity nucleotides in the flanks and a central gap of DNAs. A gapmer targeting apolipoprotein B (APOB), with five 2′-O-methoxyethyls in the flanks and a central stretch of 10 DNAs, was recently approved for the treatment of homozygous familial hypercholesterolemia,[10] and many others are in clinical development.[4] Improving on the first-generation unmodified DNA oligonucleotides, it has been generally observed that high-affinity modifications improve the potency of the antisense drug candidate, with LNA and 2′-O-methoxyethyl resulting in better potencies than 2′-O-methyl and unmodified bases.[4,11] With such an apparent proportionality between the affinity and potency of antisense oligonucleotides, it would be expected that longer oligonucleotides tend to have higher potency than shorter ones, since more nucleotides increase affinity by allowing more hydrogen bonds and additional base stackings. But this simple expectation is contradicted by experimental observations. To date, potency has been evaluated and reported for LNA-modified oligonucleotides between 12 and 20nt in length with between 2 and 5 LNAs in the flanks. When stratifying by length it has been observed that shorter oligonucleotides targeting the same target site has increased potency compared with longer versions.[7,8,11,12,13] So far, no mechanism for this seemingly counterintuitive increase in potency with decreases in affinity has been demonstrated. Suggested explanations include variations in gapsize,[7] less tendency to self-complementarity or improved pharmacokinetics of shorter oligonucleotides.[12]

In this paper, we ask if simple enzyme kinetics can account for these observations. Kinetic models have previously been successfully applied to explain and understand important actions of small RNAs or oligonucleotide-induced degradation of messenger RNAs (mRNAs). In the work by Herschlag,[14] the specificity of ribozymes and RNase H-recruiting oligonucleotides is investigated as a function of the length of their recognition sequence. As the sequence is made longer, the dissociation of mismatched duplexes is expected to become slower, and so the rate of cleavage for these mismatched RNAs will be increased relative to that of the intended target. This suggests the existence of an optimal length to achieve maximal discrimination between intended and unintended, mismatched, targets. Larsson et al.[15] examined differences in the maximal achievable knockdown of mRNAs by small interfering RNA (siRNA) and microRNA. They demonstrate that mRNAs with high turnover rates are less efficaciously downregulated since the relative contribution of an siRNA/microRNA to the degradation rate of an already rapidly recycled mRNA is small. Finally, Roth,[16] when varying parameters of a kinetic model covering uptake, trafficking, and binding of oligonucleotides, noted that cellular uptake is less important for antisense effectiveness than intracellular trafficking, oligonucleotide-mRNA binding rate, and nuclease degradation. None of these models explicitly considered the binding of RNase H to the complex of target and oligonucleotide, and thereby missed the emergence of a target- and cell type-specific optimal binding affinity, which is our primary focus here.

Results

A kinetic model of oligonucleotide-mediated RNase H recruitment and degradation of mRNA

We aim to construct a simple quantitative model of oligonucleotide action on a target RNAs. The model should be capable of addressing the observed relationship between affinity and potency. Since the experimental data available to us are from cell-based experiments, we do not consider pharmacokinetical properties such as absorption, distribution, metabolism, and excretion. Unavoidable for this class of models, some parameters will be poorly known and dynamics will be simplified and effectively renormalized.[17] If the model captures the salient features of the observations, however, it is a good indication that the mechanisms included will be of importance also in the more complex situation in vivo.

We model antisense silencing by RNase H-recruiting oligonucleotides as a four-step process (). First, the oligonucleotide (O) hybridizes to the RNA target (T). Second, the RNase H enzyme (E) binds to the complex (OT). Third, the enzyme cleaves the target to yield a complex of oligonucleotide, cleaved target, and enzyme (OCE). And fourth, the OCE complex dissociates, releasing oligonucleotide and enzyme for a new cycle, and exposing the cleaved fragments to rapid degradation by exonucleases. This system is modeled as a set of ordinary differential equations (ODEs, see Methods).

Oligonucleotides developed for therapeutics often have a residence half-life in the tissue of interest on the order of weeks.[18] Since we only model kinetic events on the order of hours, for simplicity, we assume that the total concentration of oligonucleotide, Ot, is constant. Similarly, the total amount of enzyme, Et, is assumed constant. The rate constants and initial concentrations of oligonucleotide and enzyme are set to values obtained from the literature ().

Time-resolved simulations of the model

When introducing a concentration of Ot = 0.1 nmol/l () or Ot = 100 nmol/l () at time t = 0 minutes, and simulating the model over time by numerical integration of the rate equations, we observe the expected behavior: the OT complex is formed first, closely followed by OTE. As RNA targets are degraded by the enzyme during repeated cycles, the concentrations of all compounds eventually reach a steady state, where natural as well as RNase H-mediated degradation of target is balanced by its constant production. As expected, a higher concentration of added oligonucleotide () results in larger reduction of target than a lower concentration (). Steady state is reached after ~1–2 hours as in . Experimentally, target levels are often measured after several hours or even days. Therefore, we assume steady state in the following.

Model-simulated dose–response curve

In practice, after addition of oligonucleotide to cells or tissue, target silencing is quantified by measurement of the target concentration and compared with a control situation without the oligonucleotide. In untreated cells, the steady-state concentration of target is . Since a transcript is translatable even when oligonucleotides are bound to it,[19] the total target concentration [T] + [OT] + [OTE] is the relevant quantity for evaluation of functional knockdown. This is also what is measured by quantitative reverse transcriptase–polymerase chain reaction (PCR). We use the model to calculate the steady-state ratio (Trel) of target levels in treated versus untreated (or control treated) cells as a function of Ot. This yields a dose–response curve (), with an efficacy (the maximum decrease in Trel) and potency, dependent on the choice of model parameter values (Supplementary Figures S1 and S2). Potency is defined as the inverse of the half-maximal effect concentration, EC50, which in mathematical terms is defined as:

with . According to the model, the efficacy of an oligonucleotide can be increased, by increasing the rate by which the OT complex recruits RNase H and the rate by which RNase H cleaves the target (Supplementary Figure S2). However, changing the formation rate of the OT complex has no effect on the efficacy. In general, a high concentration of RNase H improves both the efficacy and the potency, and a low target production and a low basal (non-oligonucleotide mediated) target degradation rate also improves efficacy and potency. Therefore, cells with high concentrations of RNase H and targets with low production and degradation rates are preferred for antisense oligonucleotides therapeutics.

The kinetic model predicts the existence of an optimal binding affinity

Examining EC50 as a function of the binding affinity between the oligonucleotide and the target, we see that EC50 decreases with increasing affinity to a certain point where the effect of increased affinity reaches a plateau (, dashed line). Here, the affinity is defined as the inverse of the dissociation constant between the oligonucleotide and the target (KdOT), which is in turn defined by the ratio between the dissociation and association rate of OT, i.e., . In principle, KdOT can change if either of the association or dissociation rates are changed. The association rate can be interpreted as the frequency by which a collision between the oligonucleotide and the target site will be aligned sterically in such a way as to allow formation of base pairs and thus binding. For oligonucleotides, this collision frequency is primarily limited by the accessibility of the RNA, i.e., the secondary structure of the targeted RNA.[20,21,22,23] The dissociation rate is interpreted as how often the energy of the thermal vibrations will be able to overcome the free energy of hybridization and allow dissociation. In practice, affinity is altered by adjusting the number of base pairs (by shortening or lengthening the oligonucleotide) or by adjusting the strength of each base pair (by adding or removing high-affinity modifications). Hence, we model affinity changes by modulating the dissociation rate, , while keeping the association rate, , constant. The experimental range of KdOT-values is taken from the literature.[16]

Next, we reasoned that there is an important coupling between the rate at which an oligonucleotide dissociates from an uncleaved target and the rate at which it leaves the cleaved target. After the enzyme hydrolyzes the phosphodiester backbone of the RNA target, the dissociation becomes a two-step process: either the left or right end leaves first (). The slowest of these processes will be rate-limiting for the recycling of free oligonucleotide (denoted as ). If the dissociation rate from the intact target () is decreased (e.g., due to increased affinity), must also be decreased. Moreover, since after cleavage, there will be fewer base pairs (and thus fewer hydrogen bonds and stacking interactions) to break than when dissociating from an intact target, . This is included in the model by the simple, linear relation , where 0 < α< 1. A high value of α models a case where the oligonucleotide leaves the cleaved target at a rate close to the rate at which it dissociates from the intact target. Intriguingly, irrespective of the actual value of α, the introduction of this coupling between and reveals a relationship between EC50 and affinity (, solid line), where there is an optimal affinity beyond which additional added affinity is detrimental to potency. Intuitively phrased, the higher the affinity toward the uncleaved target, the higher the affinity will also be to either parts of the cleaved target. And a high affinity toward the cleaved target results in a slow release of enzyme and oligonucleotide, which stalls the catalytic cycle and therefore reduces potency.

The ODEs behind the model (Eqs. 1–7) allow both deterministic, time-resolved, simulations (,) and derivation of algebraic solutions. We also explored stochastic simulations[24] of the ODEs. This reproduced the existence of an optimal affinity (Supplementary Figures S3 and S4), supporting the deterministic solution to the ODEs.

Experimental support for an optimal binding affinity

To investigate if the model-predicted existence of an optimal affinity could be confirmed by experiments, we reanalyzed data from two published studies evaluating the ability of different oligonucleotides to reduce luciferase (LUC) mRNA in HeLa cells[7] and glucocorticoid receptor (GR) mRNA in Hep3B cells,[13] respectively (Supplementary Material S1 for oligonucleotide data and Supplementary Material S2 for mRNA sequences). Using a validated nearest-neighbor model for LNA-modified oligonucleotide-binding to RNA (see Methods), we calculated standard free energies of binding at 37 °C, ΔG°, for those oligonucleotides in the two studies where both target knockdown efficiency and temperatures of melting, Tm (needed for the validation), were available. Plotting ΔG° values calculated in this manner against the realized target knockdown in each study (,) confirms the model's prediction of a region of optimal affinity (). We note that the model predictions of EC50 as a function of KdOT closely resembles a parabola over the experimentally observed affinity range ( and Supplementary Figure S5). To estimate an optimal ΔG° from the data we therefore fit a second-order polynomial by least squares (gray lines in ,). As seen in ,, the optimal ΔG° is estimated to −23 and −19 kcal/mol, respectively. Furthermore, to confirm the presence of an affinity optimum also for a well-established antisense target of clinical importance,[10] and independently of the length of the oligonucleotides, we synthesized and evaluated a set of compounds consisting of four 12-mers, two 13-mers, and seventeen 14-mers complementary to the APOB mRNA. The 23 oligonucleotides cover a broad range of predicted affinities. The ability of all the oligonucleotides to reduce the levels of APOB mRNA in HuH-7 cells were quantified by quantitative reverse transcriptase–PCR. Similar to the historic published data (,), this experiment also showed that oligonucleotides with a calculated ΔG° in a mid-range tended to have a better effect on APOB, than oligonucleotides on either side of this range (), with an optimal ΔG° around −18 kcal/mol. Moreover, in all three studies (–), oligonucleotides of the same length in isolation also display the affinity optimum; the 16-mers in and the 14-mers in and . From this, we conclude that the observed optimum is not caused by length per se, but rather by affinity, as indeed predicted by the model.

Finally, in both the LUC and GR transcripts, oligonucleotides are targeted to three different binding regions (starting positions in the transcripts are listed in Supplementary Material S1 for all oligonucleotides). Which region a given oligonucleotide is targeted to is indicated in , by either a gray hexagon, triangle, or square behind the point representing that oligonucleotide. This allows us to investigate whether the binding region can influence the position of the affinity optimum. In , four oligonucleotides against LUC region A (hexagons) display decreased knockdown efficiency as binding affinity increases (negative slope). In contrast, four oligonucleotides against LUC target region B (triangles) display increased knockdown efficiency as binding affinity increases (positive slope). The 13 oligonucleotides against LUC target region C (squares), span a broader range of affinities and by themselves trace out the parabolic relation. Especially for ΔG° < 25 kcal/mol (), it seems that the different binding regions explain some of the variation in the data. In , both oligonucleotides targeted to binding region A and C by themselves trace out the parabolic relation, and there are no obvious systematic variation between binding regions. In , the oligonucleotides are not focused on specific regions but tiled across the transcript (Supplementary Material S1). This may explain the larger variation around the fitted curve for this data set, which is supported in the model, since differences in accessibility of the target region can be captured as differences in , which will affect KdOT. To confirm that a single target region is sufficient to establish the parabolic relation, we designed four oligonucleotides against the same binding region in PCSK9 and evaluated reduction in transcript levels at two doses by quantitative reverse transcriptase–PCR. As seen in , these four oligonucleotides by themselves demonstrate the optimal affinity, especially for the low dose, even though the span of affinities, −17 to −21 kcal/mol, is relatively narrow (compare with –).

The optimal binding affinity is dependent on turnover rates and RNase H levels

Although the experimental data are relatively sparse, it appears as if the optimal affinity may be different between studies (–). To investigate this further, we examined the model with regard to factors, which are likely to be different between different targets or different cellular contexts (). A higher total RNase H concentration (Et) decreases the optimal affinity and increases the potency (). On the other hand, a lower Et value increases the optimal affinity and decreases the potency. A fast turnover of the target, i.e., high rates of νprod and , decreases the potency (,). However, a slow target production increases the optimal affinity, whereas a slow basal target degradation decreases it. This is similar to the findings of microRNAs and siRNAs, where it has been shown that mRNAs with a high turnover rate are more resistant to micro RNA and siRNA silencing.[15] In summary, the model predicts that the optimal affinity is influenced by the target and the context in which the reaction takes place.

Discussion

In vivo a multiplicity of factors influences the potency of an oligonucleotide. These include stability, protein binding, excretion, metabolism, uptake and diffusion into tissue, and uptake into cells and trafficking inside cells. The simple kinetic model presented here, however, provides a sufficient explanation for the experimental observation of increased potency with decreased length.[7,13] Our study suggests that for a given target (with its associated basal turnover rates) and a given cell type (with its RNase H concentration), there is an optimal affinity. With the important coupling between the dissociation rates of the oligonucleotide from intact and cleaved target ( and ), the model predicts, and the experimental data confirm that oligonucleotides deviating, in either direction, from the optimum tend to have lower potency. It is well known, and follows almost intuitively from basic principles (and our model) that oligonucleotides with too low affinity will have low potency. They therefore have to be dosed at very high and saturating amounts to bind an appreciable fraction of the target transcripts, and oligonucleotide recycling is not limiting (right-hand side of ). That oligonucleotides with too high affinity may also have reduced potency, and that it can be explained by limited recycling, has to our knowledge not been appreciated before (left-hand side of ). This insight has practical implications when designing oligonucleotide libraries for drug discovery.

Since the modeling shows that the position of the affinity optimum is parameter dependent (), our recommendation for oligonucleotide drug discovery is to initially design a set of oligonucleotides with affinities broadly covering the expected optimum. This allows the optimal affinity for that target and cellular context to be identified, as demonstrated in –. Subsequent sets of oligonucleotides can then be designed with binding affinities as close to this experimentally identified optimum as possible, keeping in mind that different binding regions may shift the position of the optimum to some extent (–). This may be due to differences in accessibility, which we speculate could influence KdOT, although measurements of target structure related to observed optimal affinities will be needed to confirm this relation. We routinely synthesize oligonucleotides with a range of calculated ΔG° values from −10 to −30 kcal/mol, and usually find the optimum between −16 and −22 kcal/mol.

While 2′ modifications to oligonucleotides are known to increase affinity, other modifications decrease the affinity, e.g., phosphorothioate backbone substitutions.[25] Previously, we considered the use of affinity decreasing modifications as a tradeoff between the loss of affinity and the gain in other beneficial properties (e.g., increased plasma protein binding and cellular uptake). However, our present results indicate that in the context of a highly affine oligonucleotide, using affinity lowering chemistry will be advantageous, since a conjugate or other modification could move a too affine compound back into the affinity optimum.

Although this study has focused entirely on RNase H-recruiting oligonucleotides, the reaction scheme may be similar for siRNAs recruiting RNA induced silencing complex, suggesting that a similar phenomena could be present in this system. Although it is also possible that other mechanisms (e.g., RNA induced silencing complexs preference for certain bases at certain positions) could have a more dominating effect on the potency.

The kinetic model is deliberately limited to contain only the reactions central to enzyme-mediated cleavage of RNA. This is important to avoid inflating the number of parameters. However, additional insights could possibly be gained by expanding the model. Possibilities include adding reactions for oligo-dimerization and self-folding (both of which can be reasonably predicted computationally)[26] or even investigate the effect of lower affinity unintended targets that could theoretically deplete the amount of oligonucleotide.[27] The source code for our model is downloadable under a Creative Commons License that enables others to easily explore these or other extensions.

In summary, we provide a simple theoretical explanation for the until now unexplained observation that long (and very affine) oligonucleotides tend to be less potent than shorter (and less affine) versions. We show that the effect is due to affinity rather than length and provide an enzyme kinetic model that explains the observation, without invoking additional complexities such as cellular uptake and protein binding. In our opinion, this demonstrates how examination of basic physical principles can inform one of the most fundamental processes of oligonucleotide based therapeutics: how to design the sequence and its modifications.

Materials and methods

Equations governing antisense oligonucleotide-mediated target degradation. The kinetic model consists of seven ODEs of the seven variables: free target (T), free oligonucleotide (O), free RNase H (E), complex of oligonucleotide and target (OT), complex of oligonucleotide, target, and RNase H (OTE), complex of cleaved target, oligonucleotide, and RNase H (OCE), and complex of cleaved target and oligonucleotide (OC):

Complex formation and breaking are denoted by rate constants k with subscripts. The target production rate is denoted by νprod and the degradation rate by . The parameter values are listed in .

Simulations. An R implementation of the FORTRAN ODE solver vode[28] is used to solve the system of seven ODEs of the model. This solver is a variable-coefficient method where the step-size is calculated at each time point. All calculations and simulations were performed in R (v.3.0.0: R Core Team (2013). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/) and is described in our ASOmodels package obtainable from https://github.com/lykkep/ASOmodel, which also contains all data from this manuscript.

Oligonucleotides. LNA-modified antisense oligonucleotides were synthesized with complete phosphorothioate backbones using standard phosphoramidite protocols on an Expedite 8900 synthesizer with a Multiple Oligonucleotide Synthesis System unit (ABI, Foster City, CA). The oligonucleotides were purified by reverse-phase high-performance liquid chromatography. Liquid chromatography-mass spectrometry (reverse phase and electrospray ionization-mass spectrometry) were used to verify compound identity and purity. For the full sequences see Supplementary Material S1.

In vitro activity screening of oligonucleotides. For APOB targeting oligonucleotides, human hepatocellular carcinoma, HuH-7, cells were transfected 24 hours after seeding with oligonucleotides at final concentration of 1 nmol/l, using lipofectamine 2000 (Invitrogen, Carlsbad, CA), according to the manufacturers instructions. Cells were harvested 24 hours after transfection and total RNA isolated using the RNeasy mini kit (Qiagen, Hilden, Germany). The PCSK9 targeting oligonucleotides were taken up unassisted in the prostate cancer cell line 15PC3, 24 hours after seeding, at final concentrations of 0.06 or 1.5 µmol/l. Cells were harvested after 3 days of uncubation with oligonucleotide and total RNA isolated using the RNeasy mini kit. In both screens, transcript levels were quantified using TaqMan assays (Applied Biosystems, Foster City, CA). The reverse transcription reaction was carried out with random decamers, 0.5 mg total RNA and the M-MLV RT enzyme (Ambion, Naerum, Denmark) according to protocol. Dependent on expression level, first-strand complementary DNA was subsequently diluted five times in nuclease-free water before addition to the reverse transcriptase–PCR reaction mixture. The Applied Biosystems 7500 real-time PCR instrument was used for amplification. Transcript levels were normalized to glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and presented as percentage-change relative to average levels in mock-treated controls. Knockdown data are available in Supplementary Material S1.

From melting temperatures to a nearest-neighbor model for binding affinity. The melting temperature is defined as the temperature at which half of the oligonucleotides are duplexed with target RNA, and is an often-used experimental measure of duplex stability, although not strictly proportional to binding affinity at physiological conditions.[29] Instead, the standard free energy of binding, ΔG° at 37 °C, is a more physiologically relevant representation of binding affinity, since KdOT is logarithmically proportional to ΔG°.[30] Also, calorimetric measurement of Tm by thermal denaturation is affected by oligonucleotide concentration and presence of cations in the buffer,[29] conditions which often differ between studies.[7,13] To allow comparisons between studies reporting Tm values, and in particular to be able to replace Tm with ΔG°, we therefore constructed a nearest-neighbor model with thermodynamic parameters taken from the literature.[31,32,33] To validate the model, we used it to calculate Tm for oligonucleotides from the two studies with historic data (,), where both knockdown data and experimentally measured Tm values were available. The residual standard error between experimental and calculated Tm is ±2.6 °C for the LUC targeting oligonucleotides[7] and ±3.5 °C[13] for the GR targeting oligonucleotides[13] (Supplementary Material S1). This is comparable with reported accuracies for nearest-neighbor models,[33] confirming the validity of the nearest-neighbor model-based calculation for these oligonucleotides. Importantly, the nearest-neighbor model also allows calculation of ΔG°, based on the same thermodynamic parameters.[34]

Author contributions

M.W.L. designed the laboratory experiments. L.P., P.H.H., and M.L. developed the model and wrote the paper. All authors have approved the final manuscript.

Figure S1. Dose–response curves for different values of Et, α, νprod, , KdOT, , KdOTE, , and (top,left to bottom,right). Figure S2. Efficacy and EC50 is plotted as functions of parameter values for Etot, KdOTE, νprod, α, , , , . Figure S3. The time-trace for the relative total target level when the model is simulated stochastically. Figure S4. Left: Dose–response curves for various values of (compare to Supplementary Figure S1, middle). Right: EC50 as a function of . A high value of corresponds to a low affinity. Figure S5. EC50 as a function of KdOT is fitted on a log-log scale to a parabola. Material S1. Oligonucleotide data. Material S2. mRNA sequences.

Table 1

Default values for the parameter-space of the model taken from the literature