Frequency-Domain Response Analysis for Quantitative Systems Pharmacology Models.

Drug dosing regimen can significantly impact drug effect and, thus, the success of treatments. Nevertheless, trial and error is still the most commonly used method by conventional pharmacometric approaches to optimize dosing regimen. In this tutorial, we utilize four distinct classes of quantitative systems pharmacology models to introduce frequency-domain response analysis, a method widely used in electrical and control engineering that allows the analytical optimization of drug treatment regimen from the dynamics of the model.

Optimizing dose and dosing regimen is arguably the most critical contribution of clinical pharmacologists to drug development. Despite significant advancements, rational and quantitative methodologies, inadequate dose, and regimen selection remains a main contributor to late‐stage drug development attrition and the need for postmarketing dose adjustments.1 Dosing regimen is typically defined as the schedule of doses of a therapeutic agent per unit of time. The conventional pharmacokinetic/pharmacodynamic (PK/PD) approach is to define a target efficacious concentration based on the two PD parameters maximum drug effect (Emax) and half‐maximal effective concentration (EC50; drug potency) and then a dose and dosing regimen based on the two PK properties clearance and volume of distribution. It has been suggested that even for drugs associated with more complex PK/PD models, these principles remain the same.1 However, recent advances in systems biology suggest that there may be an alternative approach to optimizing dosing regimen, which takes into account dosing frequency as an independent determinant of PD response. Mitchel et al.,2 for example, described that growth of yeast (Saccharomyces cerevisiae) slows down when challenged with 0.4 M KCl at a frequency of 0.125 min−1 but not at higher or lower frequencies. The authors were able to (retrospectively) rationalize this “Achilles heel” frequency using a mitogen‐activated protein kinase systems biology model. This suggests that frequency of dosing should be investigated further for modulating biological pathways and, therefore, pharmacological intervention.

That dosing frequency can significantly impact treatment success is also nicely exemplified in the re‐emergence of metronomic chemotherapy in which lower doses are administered more frequently.3, 4, 5, 6, 7 Multiple clinical trials in adult and pediatric patients with cancer suggest that such a treatment regimen could be an interesting alternative.8 Even though there have been some studies, especially in oncology, analyzing the impact of different dosing schedules on treatment success,9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 a mathematically sound or even analytical method is still lacking. Currently, pharmacometricians usually rely on trial‐and‐error methods by brute‐force simulations of only a short list of possible regimens dictated by clinical practice to find the optimal regimen for their drug and model rather than applying a more quantitative or even analytical approach up front.

Quantitative systems pharmacology (QSP) models are often described by ordinary differential equations in the time domain because the inputs (e.g., the plasma concentration of a drug or a schedule of drug administrations) and the outputs (e.g., the effect of a drug) vary in time. One fundamental aspect of these dynamic models is the timescale(s) they act on. These timescales can span a wide range, from drug‐receptor binding processes that happen within seconds to tumor growth over the course of years. Similarly, perturbations to such systems span different time scales as well and are influenced by dosing regimen and PKs. Engineers study how similar dynamic system responds to such perturbations at various time scales with frequency‐domain response analysis (FdRA).20 Therefore, they are not interested in the temporal evolution of an input signal but its harmonic content (i.e., frequency, amplitude, and phase) and how it changes when passed through a given system. The FdRA not only provides valuable insight into the dynamic behavior of a system but also enables the identification of the systems mathematical structure solely through observation of its responses to different inputs. This “black box” approach results in a so‐called transfer function that relates the inputs with the outputs, and allows the prediction of the response of the system to arbitrary inputs. Experimentally, this “black box” approach is carried out by perturbing the system at different time scales and measuring the output. When the probed system is linear and the input signal is sinusoidal FdRA is straightforward. By using, for example, multiple sound waves of different frequencies as input signals into a telephone landline and observing the output at the receiving end one can determine that the system acts as a low‐pass filter by only passing low input frequencies while filtering out high input frequencies. Although the generation of such sinusoidal signals is easy in an engineering context, only advances in microfluidics research has allowed system biologists to expose cells with similar oscillating input signals.2, 21, 22, 23, 24, 25 Systems pharmacology models, however, usually describe the interaction between drugs and a biological system.26 Thus, the inputs are usually drug administrations that do not allow for systematic probing of the system by applying a wide range of dosing regimen. Nevertheless, in order to predict which dosing regimen results in the desired system output, FdRA can be used.

With the emergence of personalized medicine in combination with the increasing prevalence of personal wearable computing devices we see a potential for uncommon but more optimal dosing regimen. Thus, the application of FdRA might help to find these optimal dosing schedules.

In this tutorial, we aim to introduce FdRA to a systems pharmacology audience by first deriving the mathematics behind it step‐by‐step with a simple example. Subsequently, we apply FdRA to four distinct and commonly used PD models (indirect response, autoregulation, precursor‐pool, and moderator‐mediated feedback models), and compare the frequency response of the linearized versions of these models with their original nonlinear version that was excited by a one‐compartment PK model with repetitive i.v. bolus administration.

A step‐by‐step derivation of FdRA

Every real‐valued continuous mathematical model that is described by ordinary differential equations can be subject to FdRA. The analytical derivation of the models’ frequency response follows the workflow given in Figure 1. Because a stable steady‐state is a prerequisite for FdRA,20 any given mathematical model first needs to undergo steady‐state analysis.27 For the steady‐state analysis here, we assume that the steady‐state input is constant (i.e., not time‐varying). If the model has at least one stable steady‐state, its linearity needs to be determined. A linear model can readily be rewritten in state‐space representation, whereas a model that is either nonlinear with respect to the states or the input must first be linearized around a stable steady‐state, which would be the one observed in patients or animals. From the state‐space representation of the (linearized) model, the transfer function can be calculated. In a last step, the transfer function is used to visualize the frequency response with the help of a so‐called Bode plot28 that related the amplitude ratio of inputs and outputs to their frequency. How a Bode plot of a nonlinear system can be determined numerically is shown as well in this tutorial.

Figure 1

Frequency‐domain response analysis (FdRA) workflow. The sequence of steps needed to perform FdRA for a given mathematical model.

An indirect response model as an illustrative example

To showcase the application of FdRA to QSP models, we select a simple indirect response model (Figure 2 a). It describes a compound that acts on the urinary bladder sphincter via stimulatory ‐adrenergic receptor, which was given to rats while a voiding volume as a physiological biomarker was followed (case study 3 in ref. 29). We denote the voiding volume with and the plasma concentration of the compound with . The drug function now is: wherein and represent the maximal effect and the potency of the stimulatory or inhibitory drug function, respectively. The differential equation describing the drug response with direct stimulation of the production is: where and denote production and loss rates of the response. Here, we should note that, unless otherwise specified, , , and are depending on time, but to improve readability we omit from the notation. According to the flowchart for FdRA in Figure 1, this model needs to undergo steady‐state analysis in the next step.

Figure 2

Indirect response model. (a) Structures of four model flavors (1 = stimulation of production, 2 = stimulation of loss, 3 = inhibition of production, and 4 = inhibition of loss) of the indirect response model. Blue arrows represent stimulation or inhibition. (b) Time course simulations of linearized model flavor 1. Input sinusoid (black) of two different frequencies ( and ) and response (pink) are shown. (c) Frequency response of all four linearized model flavors is given by the amplitude ratio for various frequencies of the input sinusoids. Schematically, low‐frequency gain, cutoff frequency, threshold frequency, and roll‐off are depicted. (d) Time course stimulations of model flavor 1. Plasma concentration as derived from pharmacokinetic (PK) model (black) for two drugs with different elimination rates ( and ) and the response (pink) are shown. (e) Frequency response of all four PK‐driven model flavors is given by the amplitude ratio for various elimination rates. Note that model flavor 1 (stimulation of production) results in the same frequency response as model flavor 3 (inhibition of production). Model parameters used in all simulations: , , , and

Steady‐state analysis

As mentioned earlier, we assume the steady‐state input to be constant (i.e., ). Thus, from the steady‐state is results in:

because the steady‐state is stable. Even though the model is linear in , the drug function (1) has a nonlinearity in . Thus, the next step is to linearize the model around the steady‐state.

Linearization

State‐space representations and, therefore, transfer functions can only be derived from linear and time‐invariant models (i.e., models that fulfil the conditions of additivity and homogeneity; see Supplementary Text). Equation (2) violates both of these conditions. This violation can be overcome if the nonlinearity in is linearized by:

This leads to the linearized model:

State‐space representation

The state‐space representation, which is just a way to represent linear time‐invariant differential equations, is now simply given by defining an output:

FdRA determines the input/output behavior of a system in response to sinusoidal inputs. Thus, we demonstrate the response of the linearized system (Eq. (3a)) to two sinusoidal inputs with different frequencies in Figure 2 b. In the left panel, an input with frequency leads to an output (i.e., response of the system) with a four times higher amplitude as compared to the unit amplitude of the input. In the right panel of Figure 2 b, an input with period into the same model results in lower output amplitude. Thus, low input frequencies amplify the response of the system, whereas high input frequencies attenuate them.

Transfer function

The transfer function now describes the input/output behavior of the system and can be determined from the Laplace transformation of Eq. (3a): as

The transfer function can now be used to collect the responses of the linearized system (Eq. (3a)) to inputs over a wide range of frequencies with the help of a Bode plot (Figure 2 c).

Bode plot

Traditionally, a Bode plot is used in electrical and control engineering to study, for example, the transmission behavior of feedback control systems. Here, we slightly modify the Bode plot to make it more accessible to a nonengineering audience. On the x‐axis, we plot the frequencies on a logarithmic scale. On the y‐axis, we plot the magnitude in base logarithmic scale. This magnitude is equal to the ratio of the output and the input amplitude (cf. Supplementary Text). Setting , , , and in (3), the Bode plot displays a constant amplitude ratio of for low frequencies and a linear decrease for high frequencies. Here, it should be noted that all amplitude ratios above entail an amplification of the input, whereas amplitude ratios below 1 represent an attenuation of the input amplitude. Engineers furthermore denote several characteristics of the frequency response,25 as visualized in Figure 2 c:

Low‐frequency gain: the degree to which constant or slowly varying inputs (e.g., long times between two drug administrations) are attenuated or amplified. Ingalls30 found that the low‐frequency gain describes the (local, steady‐state) sensitivity of the system to the input and that it reflects the steady‐state response to a step input (e.g., one i.v. bolus dose).

Cut‐off frequency: The frequency below which the system allows all components of the input to pass to the output. Alternatively, to put it differently, frequencies lower than the cutoff frequency does not perturb the system. It is the frequency for which the amplitude ratio is of the low‐frequency gain. The cutoff frequency indicates the fastest time scale the system can act on.

Threshold frequency: The frequency for which the response of the system switches from amplification to attenuation.

Peak frequency: The frequency for which the amplitude ratio is maximal. For frequency response shapes, as in Figure 2 c, the peak frequency is equal to the low‐frequency gain.

Roll‐off: The slope of the frequency response above the cutoff frequency once it reaches a straight line is called relative degree of the system, and is a measure of how quickly the output responds to changes in the input. The larger the relative degree, the slower is the response.

For the frequency response of the indirect response model in Figure 2 c, the low‐frequency gain is , the threshold frequency is , whereas the relative degree of the system follows to . The cutoff frequency can be calculated from Eq. (4) and to

In Figure 2 a and Figure 2 c, we furthermore include three additional model flavors (2 = inhibition of production, 3 = stimulation of loss, and 4 = inhibition of loss). Even though all these models are described by different differential equations and different drug functions, they all exhibit the same frequency response (i.e., the lines overlap in Figure 2 c).

Numerical FdRA

Because most QSP models are (highly) nonlinear, we contrast the analytical with a numerical approach. For that, we assume a simple one‐compartment PK model with i.v. bolus drug administration: wherein and are the plasma concentration and the elimination rate constant of the drug, respectively. The PK model in Eq. (5) is then combined with Eq. (2), i.e. the nonlinear representation of the indirect response model. The time courses for two different elimination rates of and are displayed in Figure 2 d. It should be noted that for all nonlinear models, drug administration occurs at four times the inverse elimination rate (i.e., if the drug is administered every hour). This is done to prevent accumulation in the plasma and so achieves pseudo steady‐state. As with the sinusoid‐driven linearized model, we observe higher response amplitude as compared to the plasma concentration amplitude for the lower elimination rate, whereas the higher elimination rate results in amplitude attenuation.

Performing these simulations for a whole range of elimination rate constants and measuring the plasma concentration and response amplitudes once the system reached steady‐state oscillations leads to the Bode plot shown in Figure 2 e. In comparison with the analytical frequency response given in Figure 2 c, the numerical frequency response leads to the same shape for the amplitude ratio. For low elimination rates, the amplitude ratio is nearly constant and decreases linearly for higher elimination rates. It can furthermore be seen that the four model flavors lead to three different frequency response behaviors. Although the frequency responses of the two interventions on the production are equal, those on the loss differ substantially. Over the whole range of elimination rates, but mainly for the smaller elimination rates, inhibiting the loss leads to the largest differences between plasma concentration and response amplitude, whereas stimulation of the loss leads to lowest differences. Whereas the low‐frequency gain of the analytical and numerical frequency responses differ substantially (except for the model flavor that describes inhibition of loss), the cutoff frequencies and relative degrees are comparable (Table 1). For the analytical frequency response, the threshold frequency is , whereas it is close to for the numerical frequency response.

Table 1

Characteristics of the frequency responses

Model		Low‐frequency gain		Peak frequency [ h^−1]		Cutoff frequency [ h^−1]		Threshold frequency [ h^−1]		Relative degree
type	Flavor	Anal.	Num.	Anal.	Num.	Anal.	Num.	Anal.	Num.	Anal.	Num.
Indirect response	Stimulate production	4	2.86	4	2.86	0.16	0.22	0.62	1.09	1	0.97
	Inhibit production	4	2.86	4	2.86	0.16	0.22	0.62	1.09	1	0.97
	Stimulate loss	4	2.22	4	2.22	0.16	0.31	0.62	0.93	1	0.97
	inhibit loss	4	3.99	4	3.99	0.16	0.13	0.62	1.25	1	0.98
Autoregulation	Stimulation and positive feedback	4	2.21	4	2.21	0.12	0.23	0.46	0.68	1	0.99
	Inhibition and positive feedback	4	4	4	4	0.12	0.1	0.46	0.96	1	0.98
	Stimulation and negative feedback	0.97	0.57	0.97	0.57	0.26	0.48	–	–	1	0.95
	Inhibition and negative feedback	0.97	0.89	0.97	0.89	0.26	0.24	–	–	1	0.95
Precursor‐pool	Stimulation	0	0	0.16	0.25	–	–	0.07, 0.38	0.10, 0.39	1	0.29
	Inhibition	0	0	0.16	0.15	–	–	0.07, 0.38	0.05, 0.85	1	0.37
Moderator‐mediated feedback	Stimulation	2	1.34	0.1	0.02	0.04	0.0098	0.63	1.09	1	0.98
	Inhibition	2	1.55	0.1	0.02	0.04	0.0096	0.63	1.2	1	0.98
Double moderator‐mediated feedback	Stimulation	2	1.34	0.06	0.01	0.42	0.0056	0.62	1.04	1	0.97
	Inhibition	2	1.55	0.06	0.01	0.42	0.0053	0.62	1.17	1	0.97

For all models, their flavors, and both analytical (Anal.) and numerical (Num.) frequency‐domain response analyses, the low‐frequency gain, peak frequency, cutoff frequency, threshold frequency, and relative degree are listed. For the numerical frequency response, the low‐frequency gain was calculated at an elimination rate of .

DISCUSSION

The frequency responses of the indirect response models (Figure 2 c,e) take the form of a low‐pass filter because for low frequencies the amplitude ratio is constant (i.e., those frequencies are passed to the output). High frequencies are attenuated. A direct effect model on the other hand results in constant frequency response over the whole range of input frequencies or elimination rates (Supplementary Figure S1). Thus, although all inputs are directly passed to the output for a direct effect model, certain input frequencies lead to an output amplitude modulation by the indirect response model.

Comparing the analytical and numerical frequency responses we conclude that the nonlinear systems amplify less frequent drug interventions less strongly than the linearized system. In all other frequency response characteristics, the differences are marginal. Another main difference between analytical and numerical FdRA is the inability of the analytical FdRA to resolve the different drug functions. The different model flavors only differ in the sign and location of the drug function , which is also the sole nonlinearity in the models. Thus, following linearization all model flavors are equal, which obviously leads to equal frequency responses.

If the response of a certain treatment is desired to show high fluctuations, FdRA would, therefore, advise using a low frequency dosing schedule. On the other hand, if the response should follow the plasma concentration as closely as possible, FdRA would suggest a drug dosing frequency of , whereas a minimal response amplitude would require a maximal dosing frequency.

In the following, we will apply FdRA to four distinct model structures commonly used as PD models,29 and compare the analytical approach with the numerical derivation of the frequency response.

Autoregulation models

The four autoregulation model flavors (Figure 3 a) describe autostimulatory and auto‐inhibitory processes with either stimulatory or inhibitory drug action on the loss. They can be described by: wherein , , and are the production rate, the feedback term, and the loss rate, respectively. The drug function is again given by Eq. (1). For positive and negative feedback, we defined the feedback term to be and , respectively. Herein, could be, for example, a dissociation constant.

Figure 3

Autoregulation model. (a) Structures of four model flavors (1 = stimulation with positive feedback, 2 = stimulation with negative feedback, 3 = inhibition with positive feedback, and 4 = inhibition with negative feedback). Blue arrows represent stimulation or inhibition. (b) Time course simulations of linearized model flavor 1. Input sinusoid (black) of two different frequencies ( and ) and response (pink) are shown. (c) Frequency response of all four linearized model flavors is given by the amplitude ratio for various frequencies of the input sinusoids. Note that model flavors 1 and 3 as well as 2 and 4 results in the same frequency response. (d) Time course stimulations of model flavor 1. Plasma concentration as derived from pharmacokinetic (PK) model (black) for two drugs with different elimination rates ( and ) and the response (pink) are shown. (e) Frequency response of all four PK‐driven model flavors is given by the amplitude ratio for various elimination rates. Model parameters used in all simulations: , , , , and

Analytical FdRA

For positive parameters and if we require the steady‐states to be positive as well, both unforced ( ) model flavors with positive feedback have as stable steady‐state, whereas those unforced model flavors with negative feedback have as stable steady‐state.

Thus, the linearized model in state‐space representation for the positive feedback flavors can be derived as: whereas the one for the two‐model flavors with negative feedback is: wherein and with defining as the output. Exemplarily, the response of Eq. (7a) with , , , , and , and to two sine waves with the frequencies and is shown in Figure 3 b. Although the smaller frequency leads to a fourfold amplification of the input signal, the larger frequency results in attenuation of the input amplitude.

The transfer functions of the linearized systems in Eq. (7a) are now: and where , , , and represent the models with stimulation and positive feedback, with inhibition and positive feedback, with stimulation and negative feedback, and with inhibition and negative feedback, respectively. The transfer functions in Eq. (8a) are shown in Figure 3 c over a range of input frequencies. Both the positive and negative feedback model flavors are constant for small frequencies and decrease linearly for large frequencies. Although the positive feedback model flavors possess a threshold frequency of , the negative feedback model flavors always attenuates the input independent of frequency. Thus, the low‐frequency gain of the negative feedback model flavors is less than one. Furthermore, the cutoff frequencies of all model flavors are similar (Table 1).

Numerical FdRA

All four flavors of the autoregulation model (Eq. (6)) are now excited with the one‐compartment i.v. bolus PK model given in Eq. (5). For two elimination rate constants of and , the time courses of the plasma concentration and the response of the model flavor with stimulatory drug action and positive feedback are shown in Figure 3 d. We observe that the smaller elimination rate results in higher response amplitude as compared to the plasma concentration amplitude. On the contrary, the higher elimination rate leads after a transient phase of to smaller response amplitude.

Over a range of elimination rates, the frequency response of all four model flavors is given in Figure 3 e. All amplitude ratios are constant for small elimination rates and decrease linearly for higher elimination rates. In contrast with the analytical frequency response of Figure 3 c, all four model flavors result in distinct curves. Nevertheless, only the model flavors with positive feedback possess a threshold frequency, whereas those with negative feedback attenuate the plasma concentration over the whole range of elimination rates. Furthermore, for both negative and positive feedback models, we observe that those with inhibitory drug action lead to a larger response amplitude relative to the stimulatory drug functions. Although there are some differences in the frequency response characteristics for the model flavors with stimulatory drug function, those of the model flavors with inhibitory drug function are comparable between the analytical and the numerical frequency response (Table 1).

Similar to the indirect response models, the frequency response of the autoregulation models takes the form of a low‐pass filter. However, the model flavors with negative feedback never cross the amplification/attenuation threshold. The analytical frequency response is able to distinguish between positive and negative feedback because the difference between the positive and negative feedback terms not only lead to distinct steady‐states but consequently also to distinct linear models. Its inability to resolve the stimulatory and inhibitory drug functions is based on the same issue already discussed for the indirect response models. Nevertheless, it is reassuring that, when compared to the numerical frequency response, the overall behavior is comparable, especially for the inhibitory drug function.

Precursor‐pool model

To describe tolerance and rebound, precursor‐pool models have long been used in pharmacology.31 Case study 16 in ref. 29, for example, describes the antilipolytic response of a group of healthy volunteers to an adenosine receptor agonist.32 In general, a precursor‐pool model (Figure 4 a) can be described by: wherein and represent the turnover and loss rate, respectively. Here we should note that multidimensional variables are denoted in bold (e.g., . The drug function is again given by Eq. (1). Stimulatory and inhibitory drug action now gives rise to two flavors of the precursor‐pool model.

Figure 4

Precursor‐pool model. (a) Structures of two model flavors (1 = stimulation, and 2 = inhibition). Blue arrows represent stimulation or inhibition. (b) Time course simulations of linearized model flavor 1. Input sinusoid (black) of three different frequencies ( , , and ) and response (pink) are shown. (c) Frequency response of both linearized model flavors is given by the amplitude ratio for various frequencies of the input sinusoids. Note that both model flavors result in the same frequency response. (d) Time course stimulations of model flavor 2. Plasma concentration as derived from pharmacokinetic (PK) model (black) for three drugs with different elimination rates ( , , and ) and the response (pink) are shown. (e) Frequency response of both PK‐driven model flavors is given by the amplitude ratio for various elimination rates. Model parameters used in all simulations: , , , and

Because Eq. (9) now is a two‐dimensional system, we need to use matrix notation for FdRA. As mentioned earlier, for the steady‐state analysis we assume that the input at steady‐state is constant ( ). Thus, from the steady‐state for both flavors is . The stability of the steady‐state can now be judged from the signs of the eigenvalues of the Jacobian matrix at steady‐state:

Given positive rate constants, both eigenvalues of have negative real parts. Thus, the steady‐state is stable. In order to linearize Eq. (9) around its stable steady‐states, the Jacobian matrix with respect to the model input needs to be calculated as well to:

With Eqs. (10) and (11), the linearized system can now be given in state‐space representation as: with and .

In Figure 4 b, we show the output time courses of the linearized system (Eq. (12)) for sinusoidal inputs of three different frequencies ( , , and ). Although we observe almost equal response amplitude as compared to the plasma concentration for the smallest frequency, the response amplitude is amplified for the intermediate frequency. For the largest frequency, however, the response amplitude is strongly attenuated.

For multidimensional systems, the transfer function definition in Eq. (4) needs to be rewritten in matrix form as:

Herein, denotes an identity matrix of the same dimensions as . Thus, the transfer functions for the two model flavors follow to: where and represent the model flavors with stimulatory and inhibitory drug action, respectively. Both model flavors now lead to the same bell‐shaped Bode plot (Figure 4 c). For small frequencies, we observe an increase, whereas large frequencies lead to a linear decrease. This means that the Bode plot peaks at a frequency of . This frequency response furthermore contains two threshold frequencies at and (Table 1). Thus, only within this frequency window does the plasma concentration amplitude get amplified by the system. Beyond that window, the response amplitude is lower than the plasma concentration amplitude.

Exciting the nonlinear precursor‐pool model with inhibitory drug action (Eq. (9)) with a one‐compartment i.v. bolus PK given in Eq. (5) for three different elimination rates results in the time courses shown in Figure 4 d. For a small elimination rate of , we can observe a response amplitude that is lower than the plasma concentration amplitude. An intermediate elimination rate of leads to higher response amplitude, as compared to the one of the plasma concentration amplitude. A large elimination rate of again leads to a smaller response amplitude hinting that also the nonlinear precursor‐pool model with drug PK input displays a bell‐shaped frequency response (cf. Figure 4 c). In addition, indeed, calculating the amplitude ratio for both model flavors over a range of elimination rates confirms the bell shape (Figure 4 e). For small elimination rates, we can observe an increase, whereas long elimination rates lead to a linear decrease in amplitude ratio. In contrast to the analytical frequency response, the numerical frequency response displays a plateau for intermediate elimination rates. Furthermore, we can observe than both model flavors result in distinct frequency responses with the model flavor with inhibitory drug action displaying input amplification for intermediate half‐lives. The model flavor with stimulatory drug action barely crosses the attenuation/amplification threshold. Although the frequency response characteristics are comparable between the analytical and numerical methods, the relative degree is significantly lower for the nonlinear models (Table 1 and comparing the curves in Figure 4 c,e).

The frequency response of the precursor‐pool model gives rise to a new bell‐like shape, the band‐pass filter. Although small and large frequencies result in an attenuation of the output, a small band of intermediate periods (i.e., passband) results in output amplification. The reason why the precursor‐pool model does not act as a low‐pass but a band‐pass filter can be derived from its mathematical structure (cf. Supplementary Text). Although the overall shape of the numerical frequency response is comparable with the analytical frequency response, it is interesting to observe that the passband of the numerical frequency response covers a much wider range of elimination rates. Within the passband the amplitude ratio of the nonlinear system seems to oscillate with ever decreasing amplitudes. We assume that the systems’ nonlinearities are responsible for this behavior.

The presence of a distinct peak frequency can have a significant impact on finding a suitable drug dosing schedule, as dosing at the peak frequency can lead to the most severe side effects or the most desirable effect with respect to all other drug dosing frequencies. Furthermore, for a limitation of response fluctuations both low and high drug dosing frequencies are feasible.

Moderator‐mediated feedback model

Case study 1829, 33 takes the form of a moderator‐mediated feedback model (Figure 5 a). Here, nicotinic acid inhibits the production of nonesterified free fatty acids given by in plasma. The production of free fatty acids stimulates the production of an endogenous modulator , which itself inhibits the turnover rate of the response. In general, such a system can be modeled as:

Figure 5

Moderator‐mediated feedback model. (a) Structures of two model flavors (1 = stimulation with negative feedback and 2 = inhibition with negative feedback). Green, red, and blue arrows represent stimulation, inhibition, and stimulation or inhibition, respectively. (b) Time course simulations of linearized model flavor 1. Input sinusoid (black) of three different frequencies ( , , and ) and response (pink) are shown. (c) Frequency response of both linearized model flavors is given by the amplitude ratio for various frequencies of the input sinusoids. Note that all model flavors result in the same frequency response. (d) Time course stimulations of model flavor 1. Plasma concentrations as derived from pharmacokinetic (PK) models (black) for three drugs three different elimination rates ( , , and ) and the response (pink) are shown. (e) Frequency response of both PK‐driven model flavors is given by the amplitude ratio for various elimination rates. Model parameters used in all simulations: , , , , and

Therein, , , , and are the turnover rate, the fractional turnover rate, the development of tolerance to the drug effect, and the drug function defined in Eq. (1). The presence of a stimulatory as well as an inhibitory drug function gives rise to two model flavors.

If we again assume the steady‐state input to be zero ( ) and require the steady‐state as well as all model parameters to be positive, the steady‐state for both model flavors is . The Jacobian matrix with respect to the model states evaluated at the steady‐state follows to:

Substituting , , , , and results in eigenvalues with negative real parts. Hence, the steady‐state is stable (for the given parameters). Thus, we linearize system (Eq. (15)) to arrive at the following linear models in state‐space representation:

The time courses of the linearized model with stimulatory drug action for sinusoidal inputs with three different frequencies , , and are displayed in Figure 5 b. Although the smallest frequency results in a slight amplification of the plasma concentration by the system, the input with an intermediate frequency is strongly amplified. The input with the largest frequency leads to a significant attenuation.

From Eqs. (13) and (17), the transfer functions of the two model flavors are:

Here, we denoted to be the transfer function of the model flavor with stimulatory drug action and represents the model with inhibitory drug action. The Bode plot for both model flavors is displayed in Figure 5 c. We observe that both model flavors give rise to equal frequency responses, which shows a constant amplitude ratio for small frequencies after which it peaks at a frequency of . For larger frequencies, the amplitude ratio declines linearly with a slope of . The cutoff frequency is . Furthermore, the plasma concentration amplitude is amplified for all frequencies larger than .

In order to compare the analytical frequency response of the linearized model with the nonlinear model, we again select a one‐compartment PK model with repeated i.v. bolus dosing defined in Eq. (5) to drive the nonlinear moderator‐mediated feedback model with stimulatory drug function. For three distinct elimination rates of , , and , the time courses of the plasma concentration and the response are shown in Figure 5 d. Similar to the time courses of the linearized model we observe that slow and intermediate elimination rates lead to an amplification of the plasma concentration amplitude, whereas fast elimination rates result in its attenuation. The smallest elimination rate displays lower response amplitude as the intermediate elimination rate suggesting a peak in the frequency response. In addition, indeed, numerically calculating the Bode plot of both model flavors (Figure 5 e) results in a plateau‐like peak with its maximum at . For both model flavors, this plateau lies in the amplification range of the amplitude ratio. For smaller elimination rates, we observe a constant amplitude ratio at low‐frequency gains of and for the two flavors. For larger elimination rates, a linear decrease with a slope of is present. The threshold frequency for both flavors is close to (Table 1). Last, the inhibitory drug function results in an overall higher amplitude ratio as compared to the stimulatory drug function.

The frequency response of the moderator‐mediate feedback model resembles a combination of a low‐pass and a band‐pass filter because it peaks for intermediate frequencies/elimination rates and converges to a constant amplitude ratio for small frequencies/elimination rates. In the numerically derived frequency response, we can observe similar oscillations of decreasing amplitudes as with the precursor‐pool model. Overall, the analytical frequency response as derived from the linearized model is a good approximation of the numerically simulated frequency response of the nonlinear model. This comparability is also confirmed by the similar frequency response characteristics (Table 1).

Double moderator‐mediated feedback model

Case study 19 in ref. 29 looks at the effect of an anonymized test compound upon the gene expression of specific target. Although using a stimulatory drug effect, the dynamics of fold mRNA induction was modeled by a moderator‐mediated feedback model with two moderator compartments and . In general, such a model can be described by: wherein is the stimulatory/inhibitory drug effect given in (1), the turnover rate, the fractional turnover rate of , and the fractional turnover rate for and . The presence of a stimulatory as well as an inhibitory drug function gives rise to two model flavors depicted in Figure 6 a.

Figure 6

Double moderator‐mediated feedback model. (a) Structures of two model flavors (1 = stimulation and 2 = inhibition). Green, red, and blue arrows represent stimulation, inhibition, and stimulation or inhibition, respectively. (b) Time course simulations of linearized model flavor 1. Input sinusoid (black) of three different frequencies ( , , and ) and response (pink) are shown. (c) Frequency response of both linearized model flavors is given by the amplitude ratio for various frequencies of the input sinusoids. Note that all model flavors result in the same frequency response. (d) Time course stimulations of model flavor 1. Plasma concentrations as derived from pharmacokinetic (PK) models (black) for three drugs with different elimination rates ( , , and ) and the response (pink) are shown. (e) Frequency response of both PK‐driven model flavors is given by the amplitude ratio for various drug elimination rates. Model parameters used in all simulations:

For the unforced system (i.e., without the presence of a drug and assuming that all parameters and steady‐states are positive), the steady‐state for both model flavors is To study the stability of this steady‐state we calculate the Jacobian matrix with respect to the model states of Eq. (19) to: and evaluate it at the steady‐state. For , , , , and , all eigenvalues have negative real parts. Thus, the found steady‐state is stable and can be used to linearize and express Eq. (19) in state‐space representation as:

In Figure 6 b, we display the time courses of the sinusoidal plasma concentration for , , and and the corresponding responses of the model. The slow frequency input leads to a slightly amplified response of the model, whereas the intermediate input frequency amplifies the plasma concentration amplitude fourfold. For the fastest input frequency, the plasma concentration amplitude is strongly attenuated.

With Eq. (13), the transfer function of this model is now calculated as:

The stimulatory drug effect is given by , whereas the transfer function of the inhibitory drug effect is . The Bode plot of both model flavors is displayed in Figure 6 c. Therein, we observe a constant amplitude ratio for small frequencies after which the cutoff frequency is given by . The amplitude ratio peaks at , after which it decreases linearly with a slope of . For frequencies larger than , the amplitude ratio switches from amplification to attenuation.

The moderator‐mediated feedback model with two moderators and a stimulatory drug function is now excited with a one‐compartment i.v. bolus PK model (Eq. (5)) to study the comparability between analytical and numerical FdRA. In Figure 6 d, the time courses for three elimination rates of , , and are shown. The slow and intermediate elimination rates result in an amplification of the plasma concentration amplitude, whereas the fastest elimination rate leads to an attenuation of the plasma concentration. The peak/plateau character of the frequency response is hinted at by the higher response amplitude of the intermediate elimination rate as compared to the slow elimination rate. The Bode plot confirms this observation (Figure 6 e). For the slow elimination rate, we can observe an almost constant amplitude ratio of and for both model flavors. After peaking at , the amplitude ratio reaches a plateau and declines linearly with a slope of . The threshold frequencies for both model flavors are close together. Last, the amplitude ratio of the model flavor with inhibitory drug effect is always larger than the one of the model flavor with stimulatory drug function.

Adding a second moderator compartment does not significantly change the overall shape of the frequency response that harbors a peak in the amplitude ratio for intermediate input frequencies/elimination rate. Even the frequency response characteristics are similar (Table 1). Overall, the analytical frequency response describes the response behavior of the nonlinear model well.

SUMMARY AND OUTLOOK

In this tutorial, we introduced a frequency‐domain response analysis to a pharmacometrics audience. We exemplarily applied FdRA to 14 distinct PDs models of four classes (indirect response, autoregulation, precursor‐pool, and moderator‐mediated feedback), which gave rise to low‐pass and band‐pass frequency responses. Other frequency response shapes, such as high‐pass and band‐stop, are possible as well (Figure 7 a). However, these would require the existence of the matrix (Figure 7b and Supplementary Text) that directly connects plasma concentration and response. We could imagine that such a model structure would be possible if the overall response to a drug would arise from combining one of the models we presented here (Figure 7 b gray box) with a direct effect model.

Figure 7

Frequency response possibilities. (a) Frequency responses of all possible filters (low‐pass, band‐pass, high‐pass, and band‐stop). (b) Block diagram of the state‐space representation of a linear time invariant system. Each block represents a matrix of the state‐space, whereas the block in the center is an integrator. The gray‐shaded part of the diagram highlights the model structures that were discussed in this tutorial.

To conclude, we showed that FdRA might be a helpful analytical method to not only gain insight into the dynamics of various systems but also suggests the appropriate drug dosing regimen. Additionally, in drug discovery, FdRA allows to optimize drug properties, such as elimination rate or half‐life in order to tailor the PK profile to the desired PD response. Furthermore, the comparison of analytical FdRA of the linearized models with numerical FdRA of the nonlinear models driven by a one‐compartment i.v. bolus PK model let us to believe that analytical FdRA, even though it uses drug administration in unrealistic sinusoidal inputs, provides a quick and easy way to determine which dosing frequencies have an attenuating or amplifying effect on the response of the system.

Thus, FdRA aids in the optimization of the drug dosing regimen for the desired purpose during drug development and for drug use optimization.

In the end, we should note that FdRA can not only apply to single‐input‐single‐output systems, as presented in this tutorial, but also to multiple‐input‐multiple‐output systems. In a pharmacological context, this not only allows the dose regimen optimization for combination treatments but also the analysis of multiple responses to a single or combination treatment.

In a supplementary R script, we exemplify how both analytical and numerical FdRA can be implemented.

Conflict of Interest

P.S. is funded by AstraZeneca. T.P., J.Y., and P.H.vdG. are employees of LAP&P, AstraZeneca, and Certara, respectively. As an Editor‐in‐Chief for CPT: Pharmacometrics & Systems Pharmacology, Piet van der Graaf was not involved in the review or decision process for this paper.

Supplementary Text

Click here for additional data file.

FdRA R example

Click here for additional data file.