The Well-being Model: A New Drug Interaction Model for Positive and Negative Effects

Drugs used in clinical practice usually generate two different phenomena: positive and negative effects. Negative effects manifest themselves in different ways, but they all lead to a reduction of a patient's well-being. Therefore, we consider each drug as having two main outcomes on the patient: a positive effect, resulting in well-being increase, and a negative effect, resulting in well-being decrease. The overall patient's well-being depends on the relative weights of these two effects. In the case of a combined administration of drugs, interactions might arise among the positive as well as among the negative individual effects of the drugs used. In the following, a model will be derived to describe the patient's well-being as a function of the amount of drug administered and, in the case of combined administration of drugs, of the interactions arising among them.

For any drug, the relation between the amount of drug administered and the entity of its effect, referred to as pharmacodynamics , is often modeled by a sigmoid function. For the purpose of this investigation, we refer to the pharmacodynamic model expressed as a function of the effect compartment concentration. A sigmoid curve can be adopted to describe both the positive and the negative drug pharmacodynamics. In this way, the overall patient's well-being (W ) can be defined by algebraically adding the two sigmoid curves, as shown in figure 1. Depending on the application, of course, different levels of negative effects may be tolerated. Therefore, a weighting factor (ω) for the negative effect is introduced. For one drug, the model can be summarized in the following three equations:

where W  is a parameter describing the well-being state of the patient; E  and E^N  are the drug's individual positive and negative effects; E^max  and E^maxN  are the maximum positive and negative effects achievable; EC  ⁵⁰and EC  ^50N  are the drug doses needed to reach 50% of the maximum positive and negative effect; γ and γ^N  are parameters modulating the inclination of the positive and negative sigmoid curves, respectively; C  is the drug concentration in the effect compartment; and ω is the weighting factor for the negative effects.

In this way, according to the literature,8the overall well-being curve (W ) can be interpreted as representing the therapeutic success as a function of the effect compartment concentration. Accordingly, the concentration range leading to a satisfactory value of the well-being function can be defined as the therapeutic range of the drug. Figure 1shows, for example, the positive and negative pharmacodynamics for two different drugs and the well-being function resulting from their algebraic addition. The chosen pharmacodynamic parameters for the two drugs are reported in table 1. These drugs will be referred to as drug A  and drug B . The concentrations C^A  and C^B  on this and all of the following graphs denote the effect compartment concentrations.

By modifying the pharmacodynamic parameters, curves with different therapeutic ranges are obtained. However, not all of the mathematically possible combinations of pharmacodynamic parameters represent an acceptable drug behavior for clinical purposes. Therefore, a restriction must be made to analyze only clinically significant combinations: When increasing the dose of a single drug, the positive effect is assumed to take place at lower concentrations than the negative one. In this way, the resulting curve simulates the well-being as a function of the administered dose, as observed in practice: For increasing doses, first the positive effect increases, increasing a patient's well-being; however, beyond a certain limit, negative effects arise, reducing patient's well-being and imposing a limit on the maximum dose administrable. The pharmacodynamic parameters of drugs A and B were chosen according to this criterion.

In the case of combined administration of two or more drugs, the potential interactions must be considered. Different approaches to drug interaction modeling have been developed. A recent method7is based on the use of response surfaces, where a combination of drugs is regarded as a new drug with its own sigmoid pharmacodynamic concentration-effect relation. This model will be referred to as the Minto model . In this article, the Minto model is extended to take into account also the interactions arising between the negative effects. For a better understanding of the proposed extension, the main aspects of the Minto model are briefly summarized. We will focus on a combination of two drugs, but the model is generally valid for combinations of two or more drugs. The combination of two drugs is defined by the ratio (θ) of their normalized concentrations:

where U^A  and U^B  are the normalized concentrations of drugs A and B, respectively; and θ is the drug ratio (θ ranging from 0 to 1 by definition). According to the Minto model approach, the positive effect of the combination can be modeled by a sigmoid pharmacodynamic relation where the parameters are functions of the drug ratio, of the individual pharmacodynamic parameters of each drug, and of other parameters describing the interaction effects:

where U  ⁵⁰(θ) represents the potency of the drug combination relative to the normalized potency of each drug by itself. The parameters E^max , U  ⁵⁰, and γ can be expressed as functions of the drug ratio θ by a generic polynomial of the form:

In this expression, k  is the order of the polynomial. Considering the high interpatient variability and the cost of every experimental data point, it may be not realistic to estimate a large number of parameters. Moreover, it becomes difficult to interpret the coefficients if more than the second power is present. Therefore, for simplicity, we consider the model restricted to the second order coefficients (k = 2), as already implemented in the simulations by Minto et al.  7For this case, a detailed derivation of the parameters as functions of the drug ratio can be found in literature.9The expressions simplify to

and considering that by definition U  ^50,A =U  ⁵⁰(θ= 1) = 1 and U  ^50,B =U  ⁵⁰(θ= 0) = 1, equation 11further simplifies to

where E  ^max,A , γ^A , U  ^50,A , and E  ^max,B , γ^B , U  ^50,B  are the pharmacodynamic parameters for the two drugs A and B, respectively, β^Emax , β^γ, and β^{U  50}represent the interaction parameters and E^max  (θ), γ(θ), and U  ⁵⁰(θ) are the corresponding pharmacodynamic parameters for the global sigmoid curve. In our extension, E  (θ) will be referred to as the global positive effect .

The major disadvantage of the Minto model is that it cannot be directly used to develop practical guidelines for optimal drug dosing, because it does not consider the constraints imposed by the negative effects. To overcome this, we propose to extend the model by applying the same approach to modeling the global negative effect resulting from the interactions between the individual negative effects of the drugs involved. Only a few equations have to be added to the model:

with

where all parameters have the same meaning as just described for the global positive effect. E^N  (θ) is referred to as the global negative effect . Considering the overall patient's well-being as a superposition of the global positive and the global negative effects, the complete action of the drug combination can be modeled analogously to the one drug case. By taking into account also the arising negative effects with a certain weighting factor (ω), the resulting well-being is defined as

The global positive and negative effects, as well as the resulting well-being as a function of the drug effect compartment concentrations, can be graphically described by means of a response surface. The response surface can be visualized in a three-dimensional plot, where the effect compartment concentrations of the two drugs administered lie on the two horizontal axes and the resulting effect lies on the vertical one. Figure 2shows, for example, the global positive effect, the global negative effect, and the resulting well-being for the simultaneous administration of drugs A and B, assuming that all interaction parameters are zero and for equal weight of positive and negative effects. In the following model analysis, the variations in the shape of the surface for different values of the interaction parameters will be studied to address in particular the relevant issues of uniqueness of the optimum and optimal drug dosing.

According to the model, the relation between the patient's well-being (W ) and the drug concentrations (C^A , C^B ) can be simulated for different values of the individual pharmacodynamic parameters and of the interaction parameters between the drugs considered.

In the following analysis of the model, the effect of both the pharmacodynamic parameters and the interaction parameters on the well-being surface will be investigated. Particular attention will be paid to the effect compartment concentrations required to reach the maximum well-being value.

To first analyze only the effect of the interaction parameters on the shape of the well-being surface, the pharmacodynamic parameters for the two drugs in the combination were selected (drugs A and B, as in table 1) and kept constant for each simulation. For fixed pharmacodynamics of the drugs, the effect of the interaction parameters on the shape of the well-being surface can be investigated by means of a sensitivity analysis. A sensitivity analysis consists of several simulations, each of them performed with a different value of the parameters of interest: By comparison of the simulation results, the effect of the parameters under study can be assessed. In this case, the result of one simulation is represented by the obtained well-being response surface. Two assumptions were made for the simulations:

• According to the Minto model, by assuming each drug to be able to reach the maximum effect if administered in a high enough dose, the effect can be constrained to 0 and 1 with no interaction across the E^max  (θ) parameter surface, meaning that the interaction parameters β^Emax  and β^EmaxN  can be neglected.7 

• Further, in simulations of the Minto model, no significant interactions for the parameter γ were found for the combinations studied.7Therefore, in a first approximation, we neglected it.

In this way, two simplifying equations can be added to the model:

Therefore, the simulations were performed for varying values of the interaction parameters β^{U  50}and β^{U  50N} . These parameters are referred to as β and β^N , respectively.

In the sensitivity analysis, at first, only one parameter was varied at a time. For each parameter, the two opposite cases were simulated: the parameter taking a positive or a negative value. In all simulations performed, the positive and the negative effects were weighted equally (ω=1 ). However, the analysis can be easily extended to different weighting factors.

By influencing the shape of the surface, interaction parameters also influence the effect compartment concentration range required to achieve the maximum well-being. Of course, instead of considering only the maximum value of the well-being surface, a minimal level of satisfactory well-being can be defined, considering as adequate all well-being values above the chosen minimal level. From this, the corresponding set of effect compartment concentrations leading to adequate well-being values can be derived. This is the set that should be targeted by drug infusions and is referred to as optimal concentrations range . Graphically, this set is obtained by projecting on the concentrations plane the intersection of the well-being surface with a plane at the defined minimal level of adequate well-being. For the simulations, we defined as minimal level a well-being value of 0.8, but the analysis is valid for any other chosen minimal level. Figure 3shows the derivation of the optimal concentrations range for the combination of drugs A and B, in case of zero interactions. In the sensitivity analysis, the effect of the interaction parameters on the optimal concentrations range was also analyzed.

For each simulated well-being surface, the percentage of the surface representing at least an effect of 0.8 was computed. This percentage is referred to as the Area above 0.8 . To quantify the effect of the interaction parameters, the Area above 0.8  corresponding to the simulation with β= 0 and β^N = 0 can be considered as a reference value. For each other simulated surface, the corresponding Area above 0.8  is compared with the reference value, and the percent variation, referred to as ΔArea above 0.8 , is determined according to the following equation:

The ΔArea above 0.8  can be used to quantitatively compare the different surfaces.

It could be argued that the effect of the interaction parameters on the shape of the response surface as assessed by the sensitivity analysis may depend on the chosen pharmacodynamic parameters for the two drugs. Therefore, in a second stage, the sensitivity analysis was repeated for a large number of different pharmacodynamic parameters. Further, the effect of different therapeutic ranges was analyzed in detail. Drugs with a smaller and a larger therapeutic window than drugs A and B were simulated according to the pharmacodynamic parameters reported in tables 2 and 3, respectively. These drugs will be referred to as drug C  and drug D  in the case of smaller therapeutic window and drug E  and drug F  in the case of larger therapeutic window. For a fixed value of the interaction parameters, the effect of the therapeutic range was assessed by computing the ΔArea above 0.8 , taking the surface with the smaller therapeutic range as reference.

By definition, the model is derived for effect compartment concentrations, so pharmacokinetic considerations have been neglected, so far. However, to address the issue of optimal drug dosage for clinical practice, the effect compartment concentrations entering the model must be related to the administered drug amount by the pharmacokinetics. Generally, pharmacokinetic models describe the drug distribution into the body, from the site of injection to the plasma and the different organs (or compartments). An additional model consisting of only one parameter, sometimes referred to as the link model , allows for deriving the effect compartment concentration from the plasma concentration. For the sake of simplicity, in the following analysis, the parameter of the link model is included in the parameter set referred to as pharmacokinetic parameters . If one knows the pharmacokinetic parameters of a drug, the time course of plasma and effect compartment concentrations after any intravenous administration can be computed. Therefore, deriving from the response surface the target set of optimal concentrations, one can determine the optimal drug dosage and administration pathway, to keep the effect compartment concentrations in the desired range. To analyze this, an intravenous bolus dose of 1.5 mg of drug A and 1.5 mg of drug B was simulated. Of course, all derived considerations are valid for any chosen infusion profile. A three-compartment pharmacokinetic model for drug A and a two-compartment pharmacokinetic model for drug B were assumed to be known. The chosen pharmacokinetic parameters for drugs A and B are reported in table 4. Note that no interaction was assumed between the pharmacokinetic parameters of the two drugs.

The same k  ^{e  0}was assumed for positive and negative effects. It might be argued that, in general, the effect site for positive and negative effects in the body might differ, leading to different effect compartment equilibration rates. Although this does not matter at steady state, it becomes critical during transients. In abstract modeling of drug action, however, the effect compartment is a more theoretical concept introduced to account for the hysteresis in the plasma concentration-effect relation, rather than to describe a physical location in the body. In this article, we adopt this approach. The concentration C  entering the model equations is intended as the concentration of that virtual effect compartment that minimizes the hysteresis on the well-being.

The well-being surface computed for the combined administration of drugs A and B in case of zero interactions (fig. 2) can be considered as the reference surface to analyze the effect of the interaction parameters on its shape. In figure 4, the well-being reference surface is reported again (fig. 4A), together with the results of the sensitivity analysis, where the behavior of the well-being surface is studied for varying interaction parameters (figs. 4B-I). For each surface, the value of the corresponding interaction parameters (β and β^N ) and the resulting ΔArea above 0.8  are reported. In the first four cases, only one parameter was changed at a time. Figures 4B and Dshow the effect of the interaction parameter between the positive effects, whereas the effect of the interaction parameter between the negative effects is shown in figures 4C and E. Further, simulations were performed changing simultaneously β and β^N ; these results are reported in figures 4F-Ifor the four possible combinations of positive and negative values of the two interaction parameters.

The analysis of the effect of the interaction parameter between the positive effect, β, shows that for positive values of β (fig. 4B), the surface reaches a higher well-being value already at lower effect compartment concentrations, compared with the case with zero interactions (fig. 4A). Accordingly, the percentage of the well-being surface lying above 0.8 is also increased. In case of negative values of β (fig. 4D), the opposite is true. The corresponding analysis for the interaction parameter between the negative effects, β^N , leads to the opposite result: A positive value of β^N  causes a decrease in the well-being value reached at a given combination, so that the overall Area above 0.8  is decreased. Conversely, a negative value of β^N  increases the Area above 0.8 . Combining the two effects, two cases must be distinguished: the parameters having the same or the opposite effect on the shape of the surface. For parameters having the same effect on the shape of the well-being surface (figs. 4F and G), the individual effects are enhanced. So, for example, the increase in the Area above 0.8  caused by a positive β is even higher if combined with a negative β^N  (figs. 4B and G). For parameters having opposite effects on the shape of the surface (figs. 4H and I), each individual effect is attenuated. So, for example, the increase in the Area above 0.8  deriving from a negative β^N  is lower if combined with a negative β (figs. 4E and I).

So far, the effect of the interaction parameters was derived only for a specific drug combination. The analysis was repeated for different combinations of pharmacodynamic parameters for the drugs A and B. The different pharmacodynamics affect the basic shape of the reference well-being surface, but the effect of the interaction parameters remains unchanged. Figures 5 and 6report the simulation results in the case of drug with different therapeutic windows.

In all simulations performed, although the pharmacodynamic and the interaction parameters influenced the shape of the surface, the presence of more than one maximum was never observed.

The effect of the interaction parameters is also evident on the optimal concentrations range, as summarized in figure 7Afor the effect of β and in figure 7Bfor the effect of β^N . Positive values of the parameter β enlarge the optimal concentrations range, as well as negative values of β^N , whereas negative values of β and positive values of β^N  have a shrinking effect on the same range. Figure 8shows the evolution profile of the effect compartment concentrations of drugs A and B, computed according to their pharmacokinetics. In figure 9, the optimal concentration range is combined with the evolution profile of the effect compartment concentrations of drugs A and B shown in figure 8. For simplicity, only the case of varying β (fig. 5A) is reported here.

In clinical practice, the negative effects are a limiting factor in drug dosing. The main advantage of the proposed approach is that, by including the limitations imposed by the negative effects into the expression for the well-being function, the optimal dosage can be directly defined as the one generating effect compartment concentrations lying in a range such as to lead to the maximum well-being value. Thus, the model can provide a new framework to study the behavior of the well-being function and to identify the corresponding optimal drug dosing for varying pharmacodynamics and interaction parameters.

When addressing the issue of optimal drug combinations, the basic problem is the number of possible combinations, resulting in hundreds of possible combinations that cannot be analyzed in any trial. In this case, two main approaches to the problem can be distinguished. The first is to use the information gained by analyzing a group of combinations to stepwise move toward the optimum in an efficient way. With this approach, no assumption is made on how the observed effect is generated by the combination under investigation. Examples of the application of this method to clinical questions can be found in our previous studies.3–5The second approach would be to model the observed effect as a function of the drug combination and to use this knowledge to find the combination leading to the maximum of the function. The first method has a major disadvantage: There is no guarantee that the optimum is reached.6 

Generally, a function can present more than one optimum. In this case, each optimum is referred to as local optimum , whereas only the highest of all is defined as the global optimum . In the presence of more than one optimum, a stepwise optimization procedure could lead to the identification of a local optimum only, i.e. , the identification of a set of parameters that produces a good result, but not the best one. A better knowledge of the shape of the well-being surface would help to solve this problem. The theoretical work developed in this article tries to address exactly this issue. The observed effect is modeled as a function of the drug combination, and the extensive simulations performed represent a first step in better understanding the behavior of the surface shape. One of the major results is that in all of our simulations, for the combination of two drugs, the presence of more than one maximum was never found. This would justify from a theoretical point of view the use of the direct search method for the identification of optimal drug dosing for the combination of two drugs as proposed in our previous work.4Of course, further work must be done in this direction to investigate the shape of the well-being surface also for the case of more than two drugs.

In their article, Minto et al.  7show that different interaction types (additive, synergistic, antagonistic, and others) can be simulated by the response surface approach. We agree with them that the final goal of this approach to drug interactions modeling is to characterize the response surface, more than to define how to classify the relation according to the standard definitions of synergy and antagonism. However, it would be interesting to test whether the effect of the parameters on the shape of the well-being surface as assessed by the sensitivity analysis could be consistent with an interpretation in terms of synergy and antagonism. According to Minto et al. ,7different values of the interaction parameter between the positive effects (β) can be interpreted in terms of synergistic (β > 0) or antagonistic interactions (β < 0). A detailed description of the interpretation can be found in the original article. In this extended model, for consistency, the same interpretation should be applied also to the parameter for the interaction between the negative effects (β^N ).

According to the results shown in figure 4, the effect of the parameters on the shape of the well-being surface assessed by the sensitivity analysis is consistent with the interpretation in terms of synergy and antagonism and with common clinical knowledge. Considering the interaction parameter between the positive effects, our simulations directly confirm the interpretation given by Minto et al.  7: For positive values of β (fig. 4B), the surface reaches a higher well-being value already at lower doses, consistent with the common concept of synergy. In case of negative β values (fig. 4D), the opposite is true (the surface reaches lower well-being values for a given combination of drugs), consistent with the expected effect of antagonistic interactions between the drugs.

Obviously, drug combinations leading to antagonism in the positive effect are not a relevant case in clinical practice. However, for completeness, simulations were performed also for the combination of interaction parameters describing antagonism in the positive effects. In the same way, also the interactions between the negative effects can be interpreted in terms of synergy and antagonism: A positive value of β^N  means synergy between the negative effects, whereas a negative one stands for antagonism. The effect on the shape of the well-being surface is the opposite, compared with the effect of β: Synergy between the negative effects (fig. 4C) causes them to start sooner and to reach a higher level already for lower concentrations, leading to a lower overall well-being value and to an inward curvature of the surface compared with the case with no interactions (fig. 4A). Conversely, antagonism between the negative effects (fig. 4E) attenuates the negative effects of the individual drugs, increasing the resulting well-being.

As pointed out by Berenbaum,10a therapeutic advantage may also be obtained in nonobvious circumstances, e.g. , when the positive and negative effects both show antagonism. For example, in our model, the combination of antagonism between both effects with parameters β=−2 and β^N =−2 results in an increase of the therapeutic effect (ΔArea above 0.8 =+90%), as shown in figure 4I, consistent with Berenmbaum's observation.

Figure 4Cshows that an optimal combination in a case of positive interactions between the negative effects (β^N > 0) could result in discarding one of the investigated drugs. This could explain why in a previous study, using the direct search optimization method,4clonidine disappeared from two of three optimal combinations. Namely, positive interactions between negative effects—hypotension between clonidine and bupivacaine and sedation between clonidine and fentanyl—may have resulted in a reduced well-being, leading to the elimination of clonidine during the stepwise optimization procedure.

Precise statements about the shape of the well-being surface for any specific drug combination can only be formulated if all of the required pharmacodynamic parameters for the computation of the surface are available. At this stage, only hypotheses are stated, but the test of the model on real combinations of drugs, via  previous estimation of all the required parameters, is currently under investigation. Usually, the estimation of all required parameters represents a critical point in the modeling approach with response surfaces, because of the high number of combinations to be examined in order to fit the whole surface.10An advantage of the proposed approach is that, by starting from the individual pharmacodynamics of the drugs involved, only a few interaction parameters must be estimated to completely characterize the whole response surface.

Also in the more realistic case of simultaneous interactions, the results of the sensitivity analysis are consistent with the interpretation of the parameters in terms of synergy and antagonism. Synergy in the positive effects (β > 0) tends to increase a patient's well-being. However, the entity of the well-being increase is lower if synergy in the negative effects is present as well (β^N > 0). If synergy in the negative effects is too high, the well-being increase can be counterbalanced or even changed into a well-being decrease, as illustrated by the comparison of figures 4A, B, and H, where the synergy between the positive effects first increases the percentage above 0.8 (ΔArea above 0.8 =+178%) and is then counterbalanced by the synergy in the negative effects (ΔArea above 0.8  decreasing to −54%). The same is true for antagonism in the positive effects (β < 0) combined with antagonism (β^N < 0) in the negative effects: The first one causes a decrease in patient's well-being, which can be counterbalanced by antagonism in the negative effects (figs. 4A, D, and I).

These effects remain valid for different pharmacodynamic parameters of the drugs, which only affect the basic shape of the well-being surface. Particularly interesting is the case of drugs with different therapeutic windows. The therapeutic window can be defined as the concentration range leading to a satisfactory well-being value. Thus, a drug with a larger therapeutic window will lead to satisfactory well-being values for a wider range of concentrations. When combining drugs, however, the interaction parameters influence the range of concentration leading to satisfactory well-being values. Therefore, it is interesting to analyze what is the most important factor to achieve a satisfactory well-being state for a large concentration range: strong interactions, large therapeutic window, or both? Figures 5 and 6summarize the results of this analysis. For a fixed strength of interactions, the combination of drugs with a smaller therapeutic window (figs. 5A and C) leads to satisfactory well-being values for a shorter range of concentrations than the combination of two drugs with a larger therapeutic window (figs. 5B and D). Therefore, an increased therapeutic window corresponds to an increase in the Area above 0.8 . This effect is less pronounced in the case of high interactions. On the other hand, the negative effect on the well-being due to a smaller therapeutic window can be counterbalanced by a strong desired interaction between the drugs. So, the combination of drugs with a smaller therapeutic window but a stronger desired interaction can result in a better performance than the combination of drugs with a larger therapeutic window but a smaller desired interaction, as shown in figure 6. In conclusion, according to the model, both effects are absolutely relevant in achieving a satisfactory well-being value, and no generalization can be made.

Considering the optimal combinations, the first conclusion to be drawn is that there can be more than one combination leading to the maximum well-being value. In previously published studies on drug interactions,3,4,5the application of an optimization method led to the identification of different combinations characterized by similar efficacy. Combinations of infused morphine concentration, infused ketamine concentration, and lockout interval of 1.0–1.0–8 mg/ml, 0.7–0.7–7 mg/ml, and 1.1–1.2–8 min, respectively, showed similar analgesic effect.4The presence of more than one combination with same analgesic effect indicates a plateau of the well-being surface, consistent with our simulations. Considering the less restrictive set of optimal effect compartment concentrations, as defined in figure 3, it can be stated that synergy between the positive and antagonism between the negative effects increase the optimal concentrations range, as shown in figure 7. A positive value of the parameter β enlarges the optimal area in the direction of the origin, meaning that also lower concentrations become therapeutically significant. In opposition, the parameter β^N  influences mainly the regions of higher concentrations: In case of negative values of β^N , negative effects are repressed, starting later, so that higher drug concentrations are allowed.

So far, all optimality considerations referred to optimal combinations of drug concentrations in the effect compartment. However, the more relevant question in clinical practice is how the drugs should be administered to achieve those optimal combinations of effect compartment concentrations. By combining our analysis with pharmacokinetic considerations, a first insight into the important issue of drug dosing can be developed. Considering the previously described bolus injection, the time course of the effect compartment concentrations can be computed, according to the pharmacokinetics of the two drugs involved. The plot of these concentrations on a plane having as coordinates the two effect compartment concentrations results in a loop-shaped curve, as shown in figure 8. This curve will be referred to as the effect compartment concentrations profile . Following the curve clockwise, its first part represents the increasing concentrations after the bolus dose, whereas its second part represents the decreasing concentrations according to the elimination rates of the two drugs. In our simulations, the time step between two consecutive points on the curve was 1 min. In figure 9, the effect compartment concentrations profile from figure 8is plotted together with the optimal effect compartment concentrations range derived from the well-being surface, for different values of the interaction parameter between the positive effects. Referring to the line corresponding to zero interactions (β= 0 and β^N = 0), the optimal concentrations range is achieved between 10 and 11 min. After few minutes, the concentrations reach values that are too high, leading to an inadequate well-being because of the occurrence of negative effects. After the start of the elimination phase, the concentrations will reenter the optimality zone for some time and leave it again when they become too small to provide an adequate well-being. In case of synergy between the positive effects, the optimal effect compartment concentrations range is enlarged, with the consequence that, for fixed bolus dose and pharmacokinetics, the effect compartment concentrations profile remains in the optimality region longer. For example, an adequate well-being is already achieved between 7 and 8 min after injection. The knowledge about the optimal effect compartments concentration region obtained from the pharmacodynamic model represents a basis to define optimal drug dosage and administration patterns. The optimal dosage and administration pathway can be defined as the ones keeping the resulting effect compartment concentrations inside the optimal concentrations range as long as possible. That is, the optimal administration is the one maximizing the integral of the well-being over time.

We proposed an approach to include side effects in the modeling of drug interactions for two or more drugs, allowing for a direct identification of the optimal drug combinations. Extensive simulations for the combination of two drugs showed how different drugs pharmacodynamics and interaction types influence the overall patient's well-being and the resulting optimal drug combinations range. The model is consistent with clinical knowledge and experimental results and supports from a theoretical perspective previously published approaches to the problem of identification of optimal drug combinations. The presented method improves our understanding of the characteristics of combinations used in clinical practice and provides an innovative tool to investigate drug interactions and their optimal dosage.