Biodynamics: A novel quasi-first principles theory on the fundamental mechanisms of cellular function/dysfunction and the pharmacological modulation thereof

Cellular function depends on heterogeneous dynamic intra-, inter-, and supramolecular structure-function relationships. However, the specific mechanisms by which cellular function is transduced from molecular systems, and by which cellular dysfunction arises from molecular dysfunction are poorly understood. We proposed previously that cellular function manifests as a molecular form of analog computing, in which specific time-dependent state transition fluxes within sets of molecular species (“molecular differential equations” (MDEs)) are sped and slowed in response to specific perturbations (inputs). In this work, we offer a theoretical treatment of the molecular mechanisms underlying cellular analog computing (which we refer to as “biodynamics”), focusing primarily on non-equilibrium (dynamic) intermolecular state transitions that serve as the principal means by which MDE systems are solved (the molecular equivalent of mathematical “integration”). Under these conditions, bound state occupancy is governed by kon and koff, together with the rates of binding partner buildup and decay. Achieving constant fractional occupancy over time depends on: 1) equivalence between kon and the rate of binding site buildup); 2) equivalence between koff and the rate of binding site decay; and 3) free ligand concentration relative to koff/kon (n · Kd, where n is the fold increase in binding partner concentration needed to achieve a given fractional occupancy). Failure to satisfy these conditions results in fractional occupancy well below that corresponding to n · Kd. The implications of biodynamics for cellular function/dysfunction and drug discovery are discussed.


Introduction
We proposed in our previous work [1] that time-dependent cellular function is derived from a molecular form of analog computing [2,3], in which sets of coupled ordinary differential equations and their integral solutions are modeled physically by changes in the rates of buildup and decay among the populations of biomolecular species to/from specific intra-and intermolecular PLOS  states (the hardware and software are one and the same) (Fig 1). We referred to these constructs as "molecular differential equations" (MDEs) [1], the major forms of which include: Intramolecular: Non-covalent intermolecular: Enzyme-substrate: where dγ k (t)/dt represents the rates of non-equilibrium occupancy buildup and decay of state k from one or more predecessor to one or more successor states. State transition rates are governed by adjustable barriers originating from intra-or intermolecular interactions (which we referred to as "intrinsic rates") and time-dependent changes in the concentrations or number densities of the participating species (which we referred to as "extrinsic rates") [1]. Molecular populations "flow" over time in a transient (Markovian) fashion from one specific structural state to another (Fig 1) in response to production/degradation or translocation-driven changes in the levels of the participating species.
MDEs are "solved" by cells in aggregate for the corresponding time-dependent state occupancies of the participating species {dγ k (t)/dt, k = 1, n} (e.g. dynamic ion channel states and currents, dynamic enzyme activation states, etc.), together with higher order (convergent) properties that we refer to as Γ a (t) (Fig 2). {γ k (t),k = 1, n} and Γ a (t) manifest as transient stimulus-response-driven "action potential-like" waveforms, in which the intrinsic or extrinsic rates of one or more MDEs are sped or slowed in response to Γ a (t) [1].
Our theory encompasses the major mechanisms by which: 1. MDEs are generated, powered, solved (integrated), and transduced into function at the micro-and macro-cellular levels.
2. Cellular dysfunction results from molecular dysfunction and alterations in the corresponding MDEs and Γ a (t).
3. Cellular dysfunction can be mitigated by exogenously applied MDEs, consisting of drugtarget occupancy.
In our previous work, we described the putative mechanisms by which free energy is transduced into non-equilibrium conditions and kinetic barriers (i.e. which we refer to here as "energy dynamics") [1,4,[12][13][14][15], and exemplified cellular computing by way of the cardiac AP [4] and a generic MAP kinase-phosphatase system [1]. In this work, we focus on binding as a key "integrator" of MDE systems, and in particular, the specific implications of non-equilibrium conditions on occupancy of the bound states of both endogenous molecules and drugs (i.e. the interplay between the "molecular dynamics" and "binding dynamics" branches of our theory).

Non-equilibrium conditions are both promoted and exploited by cellular systems
Rapid, transient perturbation-induced responses underlying cellular function depend on open systems, in which mass (chemical precursors, degradation products, and other substances) and energy are exchanged to/from the extracellular volume. Such responses are dampened by reversibility in closed systems exhibiting mass and energy conservation, noting that some species may undergo slow rates of buildup and decay when the overall response rate is slow. In the classical equilibrium view, covalent or non-covalent states i and j are populated according to their free energy differences ({ΔG ij }), such that the total free energy of the system is minimized. Dynamic state transitions of the participating species necessarily depend on perturbation-re-equilibration under this scenario, which is subject to the following limitations: 1. Perturbation-induced re-equilibration (i.e. equilibrium 1 ! equilibrium 2) is time-consuming, and therefore, poorly suited for rapid biological responses.
2. Optimization of complex equilibria is cumbersome and offers extremely limited evolutionary adaptability.
The resulting strain energy is proportional to the stress, analogous to that of a spring stretched away from its equilibrium state. Reversibility (e.g. A + B Ð AB) is circumvented by "tensioned" operation, in which ΔG is maintained permanently above zero via chemical or physical entry and exit of the species to/from the system over time (sources and sinks) (Fig 4). MDEs {dγ k (t)/dt, k = 1, n} are "solved" in aggregate for the corresponding time-dependent state occupancies of the participating species {dγ k (t)/dt, k = 1, n} (e.g. dynamic intra-and intermolecular states), together with higher order properties that we refer to as Γ a (t).The underlying MDEs accelerate or decelerate recursively in response to Γ a (t) according to their specific response mechanisms.

Fig 3. The six branches of biodynamics theory, which span the intra-/peri-cellular (micro-cellular) and macro-cellular levels (see text). Such processes can be
The potential energy of a system under such conditions (ΔG noneq ) is given by: where G strain is the strain energy, and G eq is the equilibrium free energy. Strained systems relax via mass action upon release of the corresponding stress, wherein ΔG noneq ! 0, and the participating species transition toward their equilibrium state distributions (i.e. [species A state i ]/[species A state j ] = K eq ). Reversibility is conveyed indirectly through concurrent (but out-of-phase) buildup and decay processes (e.g. A + B ! AB; AB + C ! A + B + C), which we refer to as "dynamic counter-balancing" ("Ying-Yang") (Figs 2 and 5). Examples include kinase-versus phosphatase-mediated phospho-transfer and inward versus outward ion channel currents. The Yin-Yang architecture serves to: 1. Prevent over-and undershooting desired γ k (t) due to the exponential behavior of molecular processes (noting that the integral solutions of MDEs are exponential functions). Yins and Yangs are necessarily maintained in an out-of-phase relationship by way of "clocking" (e.g. a phosphatase whose activation depends on a species whose buildup lags behind its target kinase).
2. Restore the initial conditions of the system.
The widespread practice of measuring biochemical (and even cellular and in vivo) effects as a function of free ligand concentration, rather than occupancy per se, reflects the extent to which the equilibrium assumption is relied upon throughout the biological sciences. However, our results suggest that time-independent equilibrium metrics of state occupancy are poor approximations to actual occupancy under cellular conditions in vivo when the rates of binding partner buildup and decay are fast.

Materials and methods
We derived two complementary analytical expressions capturing the simultaneous exponential buildup and decay of the binding site (molecular dynamics) and bound state (binding dynamics) based on a simplified second order non-competitive interaction scheme: where k on is the association rate constant, k off is the dissociation rate constant, and B free , C, and L are the free binding site, free ligand, and bound state concentrations, respectively. The rate of change of C (the mathematical equivalent of a binding MDE) is given by: Solution of this equation is facilitated by the following assumptions: 1. Molecular response is proportional solely to dynamic fractional occupancy, which in some cases, may depend additionally on the binding mode (e.g. agonist, partial agonist, inverse agonist, and antagonist in the case of receptors).
2. The total binding site level (B total (t)) is conserved instantaneously, wherein: B free ðtÞ¼ B total ðtÞ À CðtÞ ð4BÞ  (typically an enzyme), together with its substrate, product (and substrates thereof), Yang across their accessible states, as depicted by a plumbing system consisting of pipes, vessels, and valves. A given function-generating system is constituted from an interdependent collection of m such systems. Each system is maintained at non-equilibrium (analogous to "pressurization") via continuous inputs and outputs of molecular species ("fluxes"). The dynamic fluid levels within the (n � m) member vessels of the overall system, which define its instantaneous state, are solved by integration of the m MDE collections ({dγ k (t)/dt = rate of entry to state k-rate of exit from state k} i , where k = 1, 2, 3, . . ., n, and i = 1, 2, 3, . . ., m). The vessels may be connected in series or parallel. The flow rates to/from each vessel depend on the input and output valve settings (representing both the intrinsic and extrinsic free energy barriers in this simplified example), which are governed by the overall state of the system (constituting a feedback loop). (B) Dynamic counter-balancing prevents runaway exponential growth and decay of molecular species and state levels. The rates of buildup and decay consist of the net difference in Yin and Yang rates over time. The Yin-Yang balance is tipped toward the Yin during the buildup phase (left side of the blue curve) and toward the Yang during the decay phase (right side of the blue curve). The absolute rates of species or state level changes are irrelevant, except in cases where runaway behavior is functionally desirable (e.g. caspase activation). The overall buildup-decay cycle length is denoted as Λ, noting that long lived species or states likely result from slow decay rather than slow buildup (which would otherwise dampen the maximum level). Fast buildup has important implications for bound state occupancy, which are discussed later in the text.
Substituting B free (t) in Eq 4A with Eq 4B leads to: 3. That cellular systems operate in the transient regime, in which state population levels build and decay cyclically over time (although quasi-equilibrium operation is possible in some cases).
4. Buildup of B total (t) consists of the separate synthesis of the biomolecule (typically a protein) and generation of its binding-competent structural state. We further assume that state transitions occur on a timescale � the rate of synthesis of the binding site-containing species.
5. Both free and ligand-bound binding sites are lost during degradation, whereas, in practice, different rates of degradation of the bound and unbound states are possible.
6. B total (t) builds and decays exponentially via two distinct processes: a. Production and degradation of the binding site-containing species, or inter-compartmental transfer to/from the site of action, whichever is slower.
b. Formation and loss of the binding competent structural state.
7. That B total (t) builds to the equilibrium level (B 1 ) at t = 1.
8. That buildup and decay can be treated as separate sequential processes. Although B total (t) builds and decays simultaneously, net buildup occurs prior to B 1 , followed by net decay. Possible forms of the buildup term include exponential growth (Eq 6A), positive exponential decay (Eq 6B), or logistic growth (Eq 6C). Possible forms of the decay term include exponential (Eq 6D) or logistic (Eq 6E) decay (multiphasic exponential growth and decay behaviors are also conceivable): Buildup Decay where B 1 is B total at t = 1, k i and k -i are the buildup and decay rate constants, respectively, ξ is a characteristic time at which the function switches from net buildup to net decay, and t 01 and t 02 are the times at which the logistic curves reach 50% of their dynamic range during buildup and decay, respectively. We chose to approximate the buildup and decay phases of the binding site-containing species using Eq 6B (based on [22]) and 6D, and the buildup and decay phases of the binding site using Eqs 6C and 6E (Fig 6), which roughly approximates the state probability curves shown in Figs 1B and 2A. 9. k −i defines the lower limit of k off (i.e. k off is "hijacked" by the rate of binding site decay), necessitating correction of k off and K d in cases where k off < k −i : L is fixed at a constant concentration (L o ). In practice, L is expected to build and decay dynamically, which serves to further constrain the buildup of the bound state.

Scenario 1: Putative production and degradation of binding site-containing species
Substituting Eqs 6B and 6D into Eq 5 leads to: where ξ is the time required to reach equilibrium binding (conventionally referred to as the "settling time"). The net flux switches from buildup to decay at this time. ξ is defined conventionally as the time to reach 0.95 � B 1 (which we refer to as B max ), given that B 1 is only asymptotically reached over long time periods [22], whereas B max is reached over finite times.
Solving for t settling leads to: Eq 7 can be recast in dimensionless form to facilitate parameter-independent occupancy comparison, reduce the dimensionality of the system, and allow grouping of the rate terms: Solving Eq 9 for c(τ) leads to: where λ 1 = −(α + β), λ 2 = −1, and λ 3 = 0. The behavior of Eq 10A as a function of α and β can be described as follows: α + β > 0 (satisfied when α � β) ! stable conditions. α = β ! k on � L 0 = k off and cancellation of k i . Under these conditions, Eq 10A reduces to: the first term of which undergoes rapid decay at α ≳ 10 (translating to k on SS Constant fractional occupancy exists at all time points of Eq 11B, which we define as the steady-state occupancy (SSO) profile. Conversely, variable time-dependent fractional occupancy occurs at α < 10, which we define as the non-steady-state occupancy (nSSO) profile. Unlike the equilibrium case, an asymmetric relationship exists between k on and L 0 in governing dynamic occupancy, which follows from the requirement that α = β ≳ 10. The nSSO profile prevails at all non-saturating L 0 in the absence of α = β ≳ 10. For example: 1. α = β = 1 at L 0 = 1x10 -9 M, k on = 1x10 5 M -1 s -1 , k off = 1x10 -4 s -1 , and k i = 1x10 -4 s -1 (noting that only α = β is satisfied in this case).
α + β � 1 ! rapid decay of the first term in Eq 10A, which at infinite time, reduces to the equilibrium case: We assume that species or state population levels � 50% of B max make significant contributions to cellular function, translating to a "functional" buildup + decay time window: where k i and k −i likely range between ms -1 (e.g. voltage-gated ion channels) to hr -1 .

Scenario 2: putative buildup and decay of the binding competent state
Substituting Eqs 6C and 6E into Eq 5 leads to: where 2 F 1 (a,b;c;x) denotes the hypergeometric function (Fig 7). Unlike Eq 6B, Eq 6C contains a lower asymptote, which approaches zero along the negative time axis. We set B total (τ = 0) to 5% of its full maximum value, and then calculated τ 01 , τ 02 , and � x (the normalized settling time corresponding to B max = 95% of the final value), as follows: We solved Eq 14 both analytically and numerically using Wolfram Alpha (Wolfram Research, Champaign, IL, 2017) and MATLAB (Version 9.2a, The MathWorks, Inc., Natick, MA 2017), respectively. As for scenario 1, we assume that species levels or state populations � 50% of B max make significant contributions to cellular function, which leads to the following expression for Λ:

Comparison of equilibrium versus non-equilibrium occupancy
Static non-competitive, non-cooperative occupancy (γ) under equilibrium conditions can be described by the following form of the Hill equation [23]: where L o = n � K d . Fractional occupancy under the non-equilibrium SSO profile (c(τ)) is constant over time, and is also described by Eq 19. However, K d is non-equivalent under equilibrium versus non-equilibrium conditions: Non-equilibrium SSO: Achieving the SSO profile depends on in-step buildup and decay with B total (t), which in turn, depends on k on � k on SS .
Non-equilibrium nSSO:K d 0 = the upper limit of fractional occupancy at B max . nSSO occurs at all k on < k on SS , irrespective of k 0 off and K d 0 . The quasi-steady-state occupancy profile is achieved at near saturating L o levels (referred to hereinafter as qSSO). We assume the most biologically and pharmacologically desirable profile consists of SSO (i.e. where the intrinsic and extrinsic rates are tuned), which gives rise to constant fractional occupancy (B free (t)/B total (T)) over time. However, the other profiles may suffice in cases where cellular function or pharmacological/toxic effects are driven by peak, rather than sustained occupancy (e.g. ion channel blockade vis-à-vis pro-arrhythmia [4]), which necessarily coincides with C max .

Characterization of dynamic binding site occupancy under conditions of time-dependent binding site availability
We used Eqs 10 and 14 to explore the effect of decreasing Λ (and, in particular, speeding k i ) on achieving the SSO versus nSSO profile as a function of: Λ: 83,000 hr (equilibrium) to 300 ms (the approximate timescale of ion channel transitions), sampled as noted.
K d 0 (for reference): fixed at 1 nM, except where otherwise noted.

Key limitations of our approach
Our generalized analytical occupancy relationships avoid assumptions about the specific dynamic mechanisms driving binding site availability. However, such simplified relationships are subject to certain limitations compared with Markov-based simulation models tailored to specific cellular systems (e.g. as was used in our previous work [1,4]), including the potential for underestimated sensitivity of dynamic occupancy to decreasing Λ resulting from: 1) the use of fixed free ligand concentration (where ligand and binding site buildup may or may not occur in phase); and 2) neglect of competition with endogenous molecules that may also build and decay over time.

Results
Whereas static occupancy can be quantified using the Hill, Michaelis-Menten, and similar algebraic expressions, quantification of dynamic occupancy necessitates more complex timedependent differential equation-based models. However, detailed knowledge of the participating species, parameters (e.g. rate constants), and initial conditions (e.g. states and species levels) is needed to fully implement such models. We set about in this work to characterize the relationships between intrinsic and extrinsic rates and SSO versus nSSO profiles using a set of closed form dimensionless relationships that we derived for this purpose (described in Materials and methods).

Scenario 1
In this scenario, binding site lifetimes are assumed to depend on the rates of production and degradation of binding site-containing species versus formation and decay of the bindingcompetent state (i.e. transition to the binding competent state fully mirrors the buildup and decay of the species). k 0 off =k on ¼ L 0 was maintained in all cases, except as otherwise noted. We characterized the sensitivity of c(τ) to k on and k off as a function of decreasing Λ based on the following criteria: 1. The non-equilibrium threshold of Λ at which occupancy becomes kinetically-versus K ddriven (i.e. the equilibrium regime).
2. The approximate k on and k 0 off required to achieve 50% and 95% occupancy under SSO conditions (L 0 = K d 0 and 19 � K d 0 , respectively).
3. The fold-increase in L 0 required to approach the qSSO profile.
The results are summarized below and in Table 1.
The equilibrium regime. The equilibrium regime extends from Λ > 83,333 hr (k i = k −i = 1x10 -8 s -1 ), at which c(τ) is fully independent of absolute k off and k on (constrained in our study to k off /k on = K d ), converging in all cases to the SSO profile (Fig 8). Under these conditions, the bound fraction is 50% or 95% (L 0 = 1 nM or 19 nM, respectively) at all instants of time, and c (τ) exhibits the signature "saw-tooth" morphology of B total (τ). The quasi-equilibrium regime resides between 833 < Λ < 83,333 hr ( Fig 9A).
The non-equilibrium regime. c(τ) enters the non-equilibrium regime in our model at Λ ffi 833 hr (k i = k −i = 1x10 -6 s -1 ), where the SSO profile is achieved at k on � k on SS = 1x10 4 M -1 s -1 (Fig 9). Suboptimal k on results in the nSSO profile and loss of the signature saw-tooth morphology. c(τ) lags behind, and decays ahead of B total (τ), and c 50 and c 95 are no longer achieved.
Higher steady-state occupancies are achieved with decreasing Λ and k on in the range of 1x10 5 to 1x10 7 M -1 s -1 (Fig 10).
Additional simulation results as a function of Λ, k on , and k off are provided in the Supplementary Information (S1 and S4 Figs).
Achieving the SSO profile depends on fast k on , even for long-lived binding sites. We assume that long lifetimes of cognate partners are achieved via fast buildup and slow decay of the binding-competent structural state or species levels, the latter of which may consist of synthesis and degradation or inter-compartmental shuttling (noting that species dilution caused by cell growth is not necessarily well described by our exponential functions). As is apparent from Fig 11, achieving steady-state occupancy remains dependent on k on � k on SS , even at very slow binding site decay rates. At k on < k on SS , c(τ) may converge to B total (τ) far beyond the B max time point.

Scenario 2
Binding site buildup and decay in this scenario are assumed to depend on fast conformational dynamics relative to the lifetime of the binding site-containing species, which may range caseby-case from ms (e.g. voltage-gated ion channels) to hr. We assume that scenarios 1 and 2 play different roles in biodynamics (i.e. the generation of non-equilibrium conditions and MDEs, respectively). However, our overall conclusions do not depend on this assumption. As for scenario 1, we set about to characterize the effect of decreasing time-dependent binding site availability on c(τ).
The hypothetical equilibrium regime. As for scenario 1, the SSO profile is achieved in the equilibrium case in a kinetics-agnostic manner (Λ = 3,851 hr and k i = k −i = 1x10 -7 s -1 ) (S5 Fig). Under these conditions, the bound fraction is 50% or 95% (L 0 = 1 nM or 19 nM, respectively) at all instants of time, and c(τ) exhibits the same signature bi-sigmoidal morphology as that of B total (τ) (S5 Fig). The non-equilibrium regime. As for scenario 1, faster k on is needed to achieve the SSO profile as k i becomes progressively faster. The nSSO profile, and loss of the signature bi-sigmoidal c(τ) morphology, result from suboptimal k on . In such cases, c max may no longer reach Binding under non-equilibrium conditions c 50 or c 95 , and c(τ) lags behind, and decays ahead of B total (τ). c(τ) enters the non-equilibrium regime at Λ ffi 833 hr (k i = k −i = 1x10 -6 s -1 ), where a minimum k on of~1x10 5 M -1 s -1 is required to fully achieve the SSO profile (Fig 12).
Additional simulation results as a function of Λ, L 0 , k on , and k off are provided in S6 and S7 Figs.

The implications of binding dynamics for known short-lived binding sites
We next assessed the implications of binding dynamics for hERG channel blockade and LDL receptor (LDL-R) binding to wild type (wt) versus the disease-causing gain-of-function D374Y mutant form of PCSK9.
PCSK9-LDL receptor. We characterized PCSK9-LDL-R dynamic occupancy based on our scenario 1 model and the pH-dependent binding kinetics data reported in [12] (summarized in Table 2). Rapid buildup of PCSK9 can be inferred from its observed~5 min half-life [24], which is driven by zymogen activation [25] rather than de novo expression (as may be the case for many short Λ species). Negligible PCSK9-LDL-R occupancy at neutral pH follows from the exceptionally high K d relative to circulating plasma PCSK9 (which ranges between 30-3,000 ng/ml in humans (~405 pM to 40.5 nM) [26]), suggesting that extracellular binding depends on other factors. Binding and LDL-R degradation have indeed been shown to depend on heparin sulfate proteoglycans present on the extracellular surface of hepatocytes [27]. On the other hand, significant nSSO is achieved at the upper end of the concentration range for both wt and mutant forms of PCSK9 at lysosomal pH (Fig 13), suggesting that gain-of-function mutations act either through relatively small increases in dynamic occupancy relative to wt, or through some other means (e.g. by slowing PCSK9 degradation).
hERG channel blockade by non-trappable compounds. Certain hERG blockers are trapped within closed channels (which we refer to as "trappable"), whereas others are expelled during channel closing (which we refer to as "non-trappable") [4,28]. Trappable blocker occupancy builds to K d 0 (the specific time-dependence of which depends on k on ), whereas that by non-trappable blockers builds and decays during each gating cycle. The ultrashort~350 ms lifetime of the open and inactivated channel states results in k off hijacking, in which k 0 off ¼ k off ≳ k À i 3 � 0:693 s À 1 ð Þ and k off � k −i (~2 s −1 ) for the trappable and non-trappable cases, respectively [4]. In our previous work, we characterized the dynamic occupancy of a set of hypothetical non-trappable blockers using an alternate model, consisting of the O'Hara-Rudy action potential (AP) simulator, in which we replaced the Hodgkin-Huxley hERG model with a Markov state model incorporating blocker binding [4]. Here, we used the predicted proarrhythmic n � K d 0 versus n � K d as a function of k on , k off , and k 0 off as a metric of dynamic occupancy, and qualitatively compared the results with our binding dynamics models (Fig 14 ver- Table 3 versus Table 1).

Implications of binding dynamics for drug discovery
Both drug-target and endogenous time-dependent binding partner occupancy are described by similar biodynamics principles. Namely: 1. Extrinsic rates: the rates of buildup and decay of the binding site and ligand, which may vary with conditions.
2. Intrinsic rates: k on and k off , which may also vary with conditions.
3. Fractional occupancy at each instant of time (c(τ)). It is reasonable to assume that molecular response (both efficacious and toxic) depends on c(τ) � a fractional occupancy threshold during each binding site buildup/decay cycle (in general, ranging from a small fraction to nearly 100%, case-by-case). Our results suggest the paramount importance of tuning k on and k off to the rates of binding site buildup and decay, respectively, for achieving the Goldilocks zone of efficacious target occupancy and target/off-target selectivity (Fig 15). Constant maximal fractional occupancy is maintained under the SSO profile at the lowest possible L (where c(τ) depends solely on L = n � K d ). On the other hand, fractional occupancy in the qSSO regime approaches a constant level only under saturating conditions (i.e. at n >> 1) (S3 Fig). Target/off-target selectivity may be compromised with the nSSO profile, depending on the ratio of K d(target) /K d(off−target) , keeping in mind that nSSO may exist for the target, and SSO for one or more off-targets and/or competing endogenous substrates.
Tuning kinetics to the qSSO profile. Optimization of drug k off or "residence time" (i.e. t 1/2 = ln(2)/k off ) does not translate to increased occupancy when k off < k i (unless k −i is slowed in the bound state). We used Eq 10 to test the effect of slowing k off under qSSO versus SSO conditions (i.e. k on < k on_ss versus k on � k on_ss ) at constant K d (Fig 16).
Lastly, we used Eq 10 to test the effect of increasing L 0 toward the qSSO profile at k on < k on_ss (Fig 17 and S2 Fig). Binding under non-equilibrium conditions The putative effects of competition and drug pharmacokinetics (PK) on drug-target occupancy. Although we have not considered the effects of endogenous ligand competition and time-dependent drug concentration (L(t)) on dynamic drug-target occupancy, these factors can only serve to further increase k on_ss . The rates of binding site and drug buildup within the target-containing compartment (which may or may not be slower than the rate of drug buildup within the central compartment) should ideally be synchronized to ensure that L is always � K d . Fractional drug-target occupancy is decreased by L e � K i /K e in the presence of endogenous ligand competition, as is apparent from the following equation [23]: where L e , K e , and K i are the free endogenous ligand concentration and K d , and inhibitory drug K d , respectively. For example, drug SSO = 50% is achieved at L 50 = K i + L e when K i = K e , and L 50 further increases with increasing L e (or synchronization with binding site buildup) and/or Table 2. Measured binding kinetics data for wt and D374Y mutant PCSK9-LDL-R [12]. Fold differences between k on and k off for D374Y versus wt are shown in parenthesis, noting that the t 1/2 of the bound complex is similar to that of wt PCSK9, but considerably longer than that of the D374Y mutant (which is likely subject to k off hijacking). The ratio of k on SS =k on was calculated based on k on SS = 1x10 8 M -1 s -1 , corresponding to k i = 1x10 -2 s -1 (predicted from Eq 13 under the assumption that k i = k −i ).
https://doi.org/10.1371/journal.pone.0202376.t002  decreasing K e . It is further apparent that qSSO is achieved at higher L than in the non-competitive case, and that competition between endogenous ligand binding via the SSO profile would be greatly favored over inhibitory drug binding via the nSSO profile (more the reason for optimizing to k on � k on_ss ).
The binding kinetics profiles of marketed drugs are consistent with non-equilibrium drug-target binding. We set about to test the relative importance of k on_ss versus k off in the human setting using measured k on and K d data curated by Dahl and Akrud for 32 marketed drugs [29] (reproduced in Table 4). We assume that marketed drugs observe steady-state or μM blocker concentrations, respectively). The overall shape of the channel state population curves qualitatively resembles our logistic model (in which k i is very fast), noting that recovery from inactivation (the decay region of the curves) is slowed by blockade due to the smaller contribution of the hERG current to the membrane potential, which in turn, alters the response of the channel population (consistent with the recursive effect depicted in Fig 2). The blocker occupancy curves reflect nSSO occupancy at k on = 1x10 5 M -1 s -1 (where k on << k on SS , as suggested from Table 1), and the peak occupancy is far below 50% at 5 μM. Nevertheless, even transient fractional occupancy approaching the~45% level can be highly pro-arrhythmic [4]. (B) Plot of B total (τ) (gray), B free (τ) (blue) and c(τ) (gold) simulated using our logistic model (scenario 2) for a hERG blocker exhibiting K d = 5 μM (k on = 1x10 5 M -1 s -1 and k off = 0.5 s -1 ) and [free blocker] = 10 μM. k i and k −i were set to 13 s -1 and 2.1 s -1 , respectively, so as to reproduce B max at τ � 50 ms and Λ � 350 ms. B total (τ) and c(τ) are qualitatively similar to the curves in A.
quasi-steady-state occupancy for targets with fast k i and k −i , or bind to targets that build and decay slowly (i.e. kinetics-agnostic occupancy). The data can be summarized as follows: 1. k on ranges between 1.4x10 1 and 9.2x10 7 M -1 s -1 , with k on � 1x10 5 M -1 s -1 occurring in~72% of the cases, sufficient for achieving the SSO profile for binding site buildup times between 1 m ≲ t 1/2 ≲ 19 h (blue text in Table 4).
3. K d < 1 nM was achieved in~44% of the cases (red text in Table 4), which is due to k off < 3.9x10 -4 s -1 in 57% of those cases (i.e. extreme potency is not driven by slow k off ).
https://doi.org/10.1371/journal.pone.0202376.g016 [30] that the drugs were optimized on the basis of K d (noting that Folmer acknowledges the greater contribution of fast k on over slow k off to drug success [30]).

Discussion
The fundamental mechanisms by which molecular properties are transduced into cellular structure-function are poorly understood. Molecules are an unruly lot, whose states are distributed randomly in the absence of transducible energy inputs to achieve maximal configurational entropy. The relevant energy inputs consist largely of: 1) H-bond losses and gains among solvating water molecules (versus other liquid solvents) underlying intra-and intermolecular free energy barriers; and 2) continual/continuous production and decay of the participating species needed to maintain non-equilibrium conditions. Furthermore, out of phase Binding under non-equilibrium conditions dynamic counter-balancing is essential for overcoming the inherent susceptibility of exponential molecular state transitions to under-and overshooting (suggesting a putative architecture of cellular function-generating systems, consisting of Yins, Yangs, and clocks). Our biodynamics theory is aimed at capturing such principles at both the atomistic molecular and non-atomistic multi-molecular systems levels (beyond conventional mass action kinetics). Multiple levels of integrated dynamic contributions underlying cellular function, dysfunction, pharmacodynamics, and drug PK are addressed by our theory, the implications of which are summarized in the following sections.

The rates of change in binding partner levels/states are sensed via binding dynamics
The bound fraction of binding or catalytic sites can be calculated under strictly equilibrium conditions using the Hill and Michaelis-Menten equations, respectively. However, binding sites undergo transient buildup and decay driven by synthesis/degradation, translocation, or conformational state transitions to/from their binding-competent states. The lifetimes of many proteins range from single-digit minutes to hours [21] (e.g. PCSK9 is on the order of 5 min [24], and tumor cell line-derived NIK ranges from~5-30 min [31]). The lifetimes of binding competent states can range from μs (e.g. mRNA) or ms (e.g. voltage-gated ion channels) to the lifetime of the molecule itself. Here, we have shown that dynamic occupancy is greatly influenced by the time window of binding site availability, where the SSO profile is only achieved when k on is on the order of k i . We assume that the association and dissociation rate constants for endogenous partners are kinetically tuned to the dynamic ranges of B total (t) and L(t). Perturbations to partner states/levels can be sensed through binding dynamics, and specifically, kinetic tuning/de-tuning of C(t) via: 1. Modulation of the maximum time-dependent binding partner levels (B total (t)) (the "extrinsic" rates [1]).
2. Slowing/speeding of k i and k −i thereby shifting the k on threshold for achieving the SSO profile (i.e. k on SS ) or the qSSO profile.
3. Modulation of partner-specific association and dissociation rate constants (k on and k off ) (the "intrinsic" rates [1]).
Sensitivity of cellular function to the bound state depends on: 1. The dynamic range in B total (t): a. Constitutively low B total (t) levels may result in undetectable levels of the bound state.
b. Constitutively high B total (t) levels may result in reduced sensitivity of the bound state to changes in concentration of the other partner(s).
a. Very fast k on (>> k on_ss ) diminishes the sensitivity of the bound state to changes in partner concentration.
b. Very slow k on (<< k on_ss ) may result in undetectable levels of the bound state at lower partner concentrations.
We assume that cellular systems operate within the binding sensitive regime, such that the buildup and decay of binding partners and the bound states thereof remain in phase over time.
Although we have assumed a constant free ligand concentration in constructing our analytical dynamic occupancy expressions, it is apparent that this condition rarely occurs in vivo (noting that buildup and decay of the binding partners may proceed at different rates, and in an in-or out-of-phase relationship). As such, variable ligand concentration introduces potentially far greater kinetic demands on dynamic occupancy, which is likely underestimated in our approach.

The implications of biodynamics for pharmacodynamics and drug discovery
Efficacious and toxic drug levels are conventionally defined in terms of the equilibrium scenario (Eq 20), which rests on the assumptions that drug-target occupancy, drug exposure within the target compartment (typically intra-or extracellular), and target level are quasitime-independent quantities, when in fact, they are not: these conditions, efficacious or toxic occupancy, which can conceivably vary from << 50% to 95+% is expected at a free drug concentration = n � K d (e.g. n = 19 yields 95% occupancy), noting that the putative safe range of trappable hERG blocker occupancy ranges between 0-3% [4]. Conversely, achieving the qSSO profile (converging to~95% occupancy) may require considerably larger n than the SSO profile. 4. In the absence of kinetic tuning, the SSO profile is unachievable via escalation of drug concentration, although the qSSO profile is conceivable at sub-toxic exposures. This constitutes a potential source of in vitro-in vivo disconnects and clinical failures due to loss of the therapeutic index. The worst-case scenario consists of a drug that is kinetically mistuned to its target, while being kinetically tuned to one or more off-targets, and that competes for the target with an endogenous kinetically tuned partner.

5.
Residence time (i.e. ln(2)/k off below the rate of binding site decay adds no benefit to occupancy (target dynamics are not considered by other workers advocating k off optimization [32,33]).

Conclusion
In our previous work [1], we outlined a first principles theoretical framework (referred to herein as biodynamics) that explains the general mechanisms of cellular analog computing on the basis of non-equilibrium biomolecular state transitions, as well as the key implications of our theory for: 1. The transformation of chemical systems into living cells capable of maintaining non-equilibrium conditions, harnessing exponential behavior, and exploiting water H-bond energy to generate state transition barriers.
2. The origin of disease-causing cellular dysfunction (Yin-Yang imbalances), and the pharmacological mitigation thereof (insertion of exogenous MDEs capable of fully or partially restoring the normal Yin-Yang balance).
The current work is focused on the interplay between two key biodynamics contributions consisting of molecular dynamics and binding dynamics. We used a simplified analytical model that qualitatively captures these contributions to explore the possible effects of nonequilibrium conditions (namely, binding site buildup and decay) on dynamic occupancy. Our results suggest that: 1. Achieving the SSO profile as a function of diminishing binding site lifetime depends on increasingly faster k on (noting that nSSO ! the qSSO profile is achievable only at high L compared with SSO).
2. Nature has likely tuned k on to achieve SSO or qSSO binding for endogenous partners (i.e. maximal binding efficiency).
Our findings challenge the conventional equilibrium binding paradigm on which much of modern drug discovery is based, as follows: 1. Knowledge of target and binding site buildup and decay rates is essential for achieving the SSO profile on a non-trial-and-error basis.
2. Binding site lifetimes ranging from single-digit hours to milliseconds require increasingly faster k on to achieve the SSO profile (i.e. kinetic tuning versus potency optimization).
3. Increased drug exposure (i.e. PK) is a poor substitute for kinetically tuned binding, in which the qSSO profile is only asymptotically approached. As such, dose escalation during clinical trials is an extremely risky proposition from the safety standpoint. , koff = 1x10-5 s-1, and kon = 1x104 M-1 s-1, in which L0 was sampled between 1 nM and 500 nM (solid lines), and Bfree(τ) (dotted lines color coded the same as c(τ)). kon < k on SS is only compensated at sufficiently high L0, which in this case, qSSO at L0 = 500 nM.