Synergistic Interaction between Selective Drugs in Cell Populations Models

The design of selective drugs and combinatorial drug treatments are two of the main focuses in modern pharmacology. In this study we use a mathematical model of chimeric ligand-receptor interaction to show that the combination of selective drugs is synergistic in nature, providing a way to gain optimal selective potential at reduced doses compared to the same drugs when applied individually. We use a cell population model of proliferating cells expressing two different amounts of a target protein to show that both selectivity and synergism are robust against variability and heritability in the cell population. The reduction in the total drug administered due to the synergistic performance of the selective drugs can potentially result in reduced toxicity and off-target interactions, providing a mechanism to improve the treatment of cell-based diseases caused by aberrant gene overexpression, such as cancer and diabetes.


Introduction
The field of modern pharmacology aims to develop novel approaches to improve disease treatment, reduce side effects, minimize costs and enhance the efficiency of targeted therapy. These major challenges require a rational design of novel drugs and improved treatment strategies. In this direction, two of the main approaches currently being pursued involve the development of selective drugs and the design of optimal drug combination therapies.
Drug selectivity can be defined as the ability of a compound to exhibit enhanced effect towards a particular cell population in preference to others. To achieve that, a drug must be designed to target specific cellular components that are differentially expressed in two cell types. In the context of diseases that involve cells overexpressing certain genes, such as oncogenes in cancer [1], this targeting potential can be used to selectively affect only cells with increased levels of the overexpressed protein. Once selectivity is achieved, the drug can be designed to either restore normal cellular function when possible, or to trigger apoptosis of the unhealthy cells without harming the healthy cellular environment. In general, selective drugs are composed of a targeting element (TE) that recognizes and binds to the target protein, and an activity element (AE) that is directed towards the selectively targeted cells. Many of these synthetic chimeric compounds have shown good in vivo performance, and several of them have been approved by the FDA or currently undergoing clinical trials [2][3][4][5][6][7][8][9][10][11][12][13][14][15].
On the other hand, drug combination therapies have shown enhanced efficiency compared to individual drug therapy in many diseases [16], including cancer [17,18] and HIV's [19]. The interaction between drugs in multicomponent therapies is a complex and multi-scale problem [20] that requires full characterization of the direct and indirect molecular aspects of the interaction, which are often unknown. Due to this, experimental studies and discoveries of successful drug combinations are often based on empirical intuition and trial-and-error approaches. In general, drug interactions can be classified based on their effect when combined, compared to their effect when applied alone. Drugs that do not interact with each other, or are mutually exclusive by competing for the same target are considered as additive [22]. This basically means that the lower concentration which produces a certain effect corresponds to the most potent drug, and there is no gain due to the combination of the two drugs. On the other hand, antagonism occurs when one of the drugs mitigates or counteracts the action of the other, i.e, the combination is always less effective than the single agents at the same concentration. Finally, synergism occurs when the combination of both drugs is more effective than each agent separately at the same total concentration, i.e., one of the agents enhances the actions of the other [21]. This can occur either via direct interaction, i.e, one drug increases the bioavailability of the other, or indirectly, i.e, the two drugs cooperate on targets on the same or different pathways involved in the same process [23]. Thus, the total concentration of drug administered to achieve a certain effect is reduced, which potentially also reduces side effects, drug resistance and undesired off-target interactions.
In the context of selective drugs, synergism and antagonism can be also defined in terms of the enhanced or reduced selective potential of the two drugs when combined [24], i.e, their ability to target selectively a specific cell population, compared to their selective potential when applied individually. In this way, two drugs are synergistic if their combination is more selective than the two drugs acting alone at the same total concentration. Here, we explore the mechanism of interaction between selective drugs in combination from a theoretical perspective. To do that, we develop a population model where two sets of cells expressing different levels of a target molecule are treated with different concentrations of two drugs simultaneously. In principle, these two drugs can be monomeric non-selective ligands (i.e., they do not differentiate between healthy and unhealthy cells), or chimeric ligands, composed of an AE and linked to a TE, allowing them to selectively target unhealthy cells, leaving the healthy environment undamaged.
Two different approaches are taken into account: first, we analyze the effect of combinations of two different chimeric ligands when applied simultaneously to an heterogeneous population of cells; next, we combine the effect of individual chimeric drugs based on the Loewe additivity model [22]. Both models predict that drug combination of selective drugs is synergistic in terms of their selective potential. Finally, we introduce phenotypic inheritance in the cell population to show that both selectivity and synergism also occur in a context where the amount of target proteins of the daughter cells depends on the mother cell, such as in diseases caused by mutations in specific genes. Our results show that the concentration to obtain a desired selectivity can be minimized by simultaneous treatment of selective drugs.

Models
To analyze the effect of selective drug combinations in a multicellular approach, we develop a mathematical framework where we allow two asynchronous populations of cells with two distinct average number of target molecules to proliferate for a given time. Cells are treated with different concentrations of monomers and chimeric drugs alone or in combination, to then monitor and compare the dynamics of proliferation of the two cell populations. The dynamics of the effect of chimeric drugs at the cellular level is calculated based on a chemical kinetics model for the ligand-receptor interaction [25]. The model used is an extension of our previous contribution to the study of the dynamics of chimeric ligands, where we develop a mathematical framework to predict the selective potential of chimeric drugs, based on the affinity of both AE and TE subunits of the ligand towards their targets (activity element receptor, AER, and targeting element receptor, TER, respectively), the concentration of the target molecules and the linker length between AE and TE in the chimera [25]. The model is rewritten to take into account the simultaneous interaction of two chimeric ligands, resulting in the following set of interactions: where R i corresponds to the i receptor (i = 1,2) and L j corresponds to each of the two different ligands (j = 1,2) used in the combined treatment. Each ligand L j is composed of a AE and TE, and it can bind to R 1 or R 2 via reaction 1 to give an intermediate complex C i,j (Fig. 1A-C). These intermediate complexes facilitate reactions 2 and 3 by originating a local concentration of the free subunit of the ligand L j in the vicinity of the complementary receptor, to generate the complex C 3,j (Fig. 1D). The coupling (k c i;j ) rate constants in reaction 2 and 3 are calculated as follows: where the diffusive rate constant k diff i is modulated by the diffusion D i of the receptors R i at the membrane as [26,27]: being A the average cell surface area, b i ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A=ðpR i Þ p corresponds to half the average distance between R 1 and R 2 , and a is the linker length. Effective affinity and dissociation rates for the reactions that take place at the membrane are calculated as [26]: where N av is Avogadro's number and V i = A Á (h i + a) is the effective reaction volume for the second binding event, assumed as a spherical gasket above the cell surface where the free subunit gets distributed after the first binding event (see [25]), being h i the height of R i above the cell surface. g i;j ¼ k on0 i;j =ðk diff i þ k on0 i;j Þ corresponds to the capture probability factor for receptor R i and ligand L j , explained in detail in [26].
For a given constant concentration of both ligands L j , the equations are solved for each individual cell in the population, based on its amount of R 1 and R 2 receptors, identified as TER and AER, respectively. The maximum value of AER-AE complexes formed in each cell is then correlated with the physiological response produced by the AE using experimental dose-response curves (this correlation is a multi-step process explained in detail in Ref. [25]). Typical doseresponse curves are often fitted to a four-parameter sigmoidal [25], such as: where the physiological response R(L) for a given drug concentration L is characterized by its maximal D and minimal A asymptotes, B is the slope parameter of the curve, and EC 50 is the half-maximal effective concentration of the ligand, that is, the inflection point of the curve. S1A Fig. shows a schematic representation of the workflow used to solve the model equations and obtain the dynamics of growth of the heterogeneous cell population. As a numerical solution, the model is informed with data from a synthetically designed chimeric ligand composed of the Epidermal Growth Factor (EGF) as TE and different mutants of Interferon alfa-2a (IFNα-2a) as AE [15]. Thus, the apoptotic effect triggered by IFNα-2a stimulation is directed towards cells overexpressing the Epidermal Growth Factor Receptor (EGFR). The physiological response of the cells to the treatment corresponds to the apoptotic effect induced by IFNα-2a, measured experimentally as the percentage of surviving cells after 60 hours of treatment.
Since EGFR is an oncogene overexpressed in a number of tumor cells [28], this chimera can be potentially used to selectively target cancer cells without affecting the healthy surrounding tissue. Different mutants of the IFNα-2a molecule are tested as monomers (M wt , M 1 , M 2 and M 3 ), and as AE's in chimeric configuration, identified here as Ch wt for the chimera composed of the wild type version of the IFNα-2a linked to EGF, and Ch 1 , Ch 2 for the experimentally available mutants of IFNα-2a with reduced affinity towards the IFNα-2a receptor linked to the targeting element EGF [15]. Other potential chimeras composed of EGF linked to IFNα-2a mutants with decreasing affinity towards the AER combined with EGF, named Ch 3 and Ch 4 , are included in the analysis ( Table 1 shows the dissociation constants for each IFN monomer).
To mimic the experimental conditions, cells in the population are allowed to proliferate for 60 hours in the presence of the drug treatment, and the physiological response of each cell to the treatment depends on the amount and efficiency of each ligand, the amount of TER and AER receptors expressed, and the exposure time to treatment. Given that the physiological response of the cells is apoptosis, we set the decision between survival or death for each cell in the population as follows: for a given physiological response (0 < R(L) < 1), the probability of undergoing apoptosis at every time point is computed as θ = (1 − R(L)) Δt/T , being T the total length of the experiment and Δt the time step in the simulation. Then, a random number (0 < γ < 1) from a uniform distribution is assigned for each cell in the population and compared to the value of θ at every time step. If γ ! θ, the cell survives. On the contrary, if γ < θ, then the cell dies and it is no longer considered in the simulation. A cell division occurs when the age of a given cell reaches the numerical value for the cell cycle length assigned to that particular cell. This value is obtained from a gamma distribution with mean m = 27.5 hours [29,30] and standard deviation of 2 hours. The amount of surface receptors are also gamma distributed [31], with a mean value obtained from experimental data [15] (see Table 1) and a coefficient of variation of 0.3, to mimic cell-to-cell variability in both populations. Mean values of the final cell numbers are obtained from 10 independent runs of the model. Numerical solution of the model equations and other calculations are performed using a in-house Matlab script (code available upon request).

Selectivity of chimeric drugs versus monomers in a cell population
The model described above is used to illustrate the effect of monomers versus chimeric ligands in a heterogeneous cell population. Fig. 2 illustrates the dynamics of growth of healthy (blue curve) and unhealthy (red curve) cell populations under nonselective monomers versus selective chimeric ligand treatment. To characterize and compare its selective potential, we define a threshold based on the amount of cells of both populations that remain after 60 hours of treatment (i.e., the duration of the experiments in [15], where the dose-response curves and other experimental data are obtained). Thus, a given treatment is considered as efficient when the number of unhealthy cells does not increase, while the population of healthy cells grows to at least 80% of its potential size. These two threshold values are marked in Fig. 2A-I as dashed red and blue horizontal lines for the unhealthy and healthy cells, respectively. These threshold values will be used to categorize the selective potential of a given treatment.
Low concentrations of M wt are harmless to both cell populations, which can grow exponentially ( Fig. 2A). Intermediate concentrations (Fig. 2B) have a much stronger effect in healthy cells than in unhealthy cells, due to higher resistance to IFNα-2a treatment in the unhealthy cell population (reported experimentally in [15], where authors hypothesized that this effect is mainly due to the anti-apoptotic potential of the EGFR overexpression [32] that may counteract the effect of IFNα-2a stimulation). Higher concentrations of the monomer are able to reduce the number of unhealthy cells in the system, but affecting the healthy population even more (Fig. 2C). This type of response is similar for all mutants of the monomer of IFNα-2a, as shown in the S2 Fig  AER expressed by each individual cell, the model predicts that a number of cells with high values of AER does not undergo apoptosis after 60 hours of treatment. This is due to the fact that the physiological response of a given cell depends on the time that it has been under treatment and, since cells are continuously being born during the simulation, at t = 60 hours some recently born cells may not have been enough time under the influence of the drug to trigger apoptosis.
The above data evidences that selectivity in terms of the threshold of 80% can only be achieved using low affinity mutants of the AE subunit of the chimera, which are only efficient at very high concentrations of ligand (around 2 nM concentration for Ch 2 ). This concentration representes a 4X increase compared to the concentration for Ch wt = 0.5nM, which is the minimum concentration required to prevent the expansion of the unhealthy cell population at the expense of affecting strongly the healthy cells. This increase is also reported experimentally in [15], where the minimum concentration of Ch 2 to prevent the unhealthy cell population to expand leaving the 80% of the healthy surrounding undamaged is 1.5nM, 3 times higher than the value of 0.5nM of Ch wt that prevents the growth of the unhealthy cell population. Unfortunately, this higher doses of drug required to achieve selectivity can result in the emergence of other potential undesired effects, such as toxicity, or increased off-target interactions [33]. Therefore, strategies to reduce the total drug concentration for a given selective effect are relevant. In the next section, we show how combinations of selective chimeric ligands can reduce the concentration of total drug administered maintaining the selective potential.

Synergistic interaction of selective chimeric drugs
Mutant monomers of the same molecule act against the same target, therefore they behave as mutually exclusive and their interaction when combined is additive by definition. In this way, the lowest concentration to obtain a given effect always corresponds to the monomer with stronger binding affinity towards its receptor. Fig. 3A-C corresponds to simultaneous treatment of a fixed concentration of M wt = 0.5nM with the minimal concentration of the mutants of IFNα-2a able to affect at least 20% of the unhealthy cells. According to their additive interaction, none of the potential combinations tested allows us to reduce the total drug concentration of IFNα-2a administered. In addition, none of the multiple combinations tested is able to mitigate the strong effect that the monomers exhibit towards the healthy cell population, as illustrated also in the next section.
Multi-drug treatment using selective drugs is shown by Fig. 3D-I, where we plot the dynamics of the two cell populations at the minimal concentration of total drug required to achieve the selectivity threshold for different combinations of the chimeras. Combinations of the poorly selective Ch wt with chimeras Ch 2 , Ch 3 and Ch 4 are able to mitigate the expansion of the unhealthy cell population while meeting the 80% survival threshold of the healthy cells (Fig. 3D-F). The threshold is also achieved when combining the rest of the chimeras with Ch 1 , as shown in Fig. 3G-I.
Bars in each panel represent the minimal concentration of total drug to achieve the threshold of selectivity for each combination of ligands, compared to the same ligands as single treatment. We observe a slight reduction in the total concentration used for the combinatorial  Table 2). Overall, the gain in performance of the dual treatment is more evident when combining a poorly selective ligand as Ch 1 with a low affinity but highly selective chimera, such as Ch 3 . In this situation, we can achieve the threshold of optimal selectivity with a 50% reduction in the total drug concentration administered, compared to the concentration of Ch 3 as individual treatment.

The effect of combinatorial treatment can be estimated based on the effect of single treatment strategies
Detailed analysis of the drug interaction for each combination of two given ligands, as performed in the previous section, requires extensive computational resources. To overcome this, we use an additional approach to calculate the effect of a combination of drugs based on their effect when applied individually, using the following equation [22]: where L c,1 and L c,2 correspond to the concentrations of the two drugs that produce a given effect when applied together, and L 1 and L 2 are the concentrations that induce the same effect when applied alone. The interaction index I depends on the type of interaction between the two drugs: synergistic (I < 1) additive (I = 1) or antagonistic (I > 1). In our particular case, since both ligands share common binding sites (i.e., they are mutants of the same molecule with reduced affinity), they behave as mutually exclusive and therefore, they follow the principle of Loewe additivity so the interaction index I in Eq. 10 is set to 1 [22] (a detailed analysis of the general equation for drug interaction in conditions of drug additivity, synergy or antagonism can be found as S1 Text). Next, the Loewe additivity model is applied to drugs that induce a typical dose-response curve [25], as described in Eq. 9. Thus, we solve Eq. 9 for L and substitute into Eq. 10 for both ligands (L = L 1 and L = L 2 ), to obtain: This equation can be solved numerically to obtain the physiological response R(L c,1 , L c,2 ) for each potential combination of L c,1 and L c,2 . The typical shape of the curve is shown in S4A Fig. Therefore, Eq. 11 allows us to calculate directly the effect of the two drugs when applied simultaneously, significantly reducing the computational cost of the process. S5 Fig. plots the dynamics of several combinations of monomers and chimeras using this method (to be compared with Fig. 3, computed using the simultaneous drug stimulation of the system to show that both methods produce equivalent results). This simplified method allows us to compute the effect of any combination between two given drugs to develop isobolograms representing the final number of cells, after 60 hours of combined treatment for each combination of monomers (Fig. 4A-B) and chimeras (Fig. 4D-E).
Finally, to compare the selective effect of multiple combination of ligands, we define the performance P(L c,1 , L c,2 ) of a given treatment as: where N Ã f ðL c;1 ; L c;2 Þ healthy and N Ã f ðL c;1 ; L c;2 Þ unhealthy correspond to the final number N f of healthy and unhealthy cells respectively, after 60 hours of combined treatment (i.e., the final point of the curves in Fig. 3) normalized to obtain values for the performance P(L c,1 , L c,2 ) between −1 (minimal selectivity of treatment, i.e., 0% survival of the healthy cells) and 1 (optimal selectivity, i.e., !80% survival of healthy and no growth in the unhealthy cell population). A workflow scheme of this approach is shown in S1B Fig.  Fig. 4A-B illustrates the isobolograms for any concentration of monomers M wt and M 1 for healthy (Fig. 4A) and unhealthy cells (Fig. 4B) for a range of concentration values. The performance map (Fig. 4C) evidences that monomer combinations are not selective (i.e, P 0 at any concentration). Isobolograms for Ch wt and Ch 1 combinations are shown for healthy (Fig. 4D) and unhealthy (Fig. 4E) cell populations. The performance map (Fig. 4F) illustrates that the combination shows regions of positive performance (P > 0), i.e., regions where the combinatorial treatment acts selectively towards the unhealthy cell population.
Performance colormaps for other combinations of chimeric drugs are shown in Fig. 5, where regions in which the threshold of selectivity is achieved are marked in dark red. Simultaneous treatment of Ch wt with Ch 2 , Ch 3 and Ch 4 show that, for each combination, optimal selectivity can be achieved at slightly lower concentrations when using the two drugs simultaneously, compared to the same drugs acting alone (i.e., [Ch wt ] = 0 in each panel). This is more evident when combining Ch 1 with the other chimeras (Fig. 5D-F), where the optimal selectivity threshold is achieved at significantly lower concentrations using two selective drugs instead of one (numerical values for the minimal concentration for the threshold of selectivity as well as the reduction in total final concentration due to the combinatorial treatment are shown in Table 2).

Synergistic interaction for selective chimeric drugs in cell populations with heritability
Previous simulations assumed that both healthy and unhealthy cells are, as a first approximation, phenotypically different, with the expression levels of TER and AER obtained from gamma distributions. Therefore, the amount of receptors expressed by a daughter cell depends on its cell type, but it is independent on the amount of receptors expressed by the mother cell. In other potential scenarios, the difference in phenotype between healthy and unhealthy cells can be caused by genetic mutations, and therefore, the amount of receptors expressed by the mother cell is inherited by the daughter cells. In these situations, a given treatment can become inefficient, and it can potentially act as selective pressure, acting more strongly over weak cells and ultimately inducing resistance to treatment in the population. This scenario has been explored extensively in vivo and in silico, and is one of the main causes of the short-lived response of targeted therapy in cancer [18]. Recently, it has been shown theoretically and experimentally that dual treatment strategies can dramatically reduce the possibility of development of resistant cells, resulting in long-term disease control compared to single drug treatment, or even sequential drug treatment [18].
To test the performance of dual selective treatment in the context of genetic inheritance of the amount of receptors, we set the average amount of TER and AER expressed by a given daughter cell as directly given by the amount of receptors expressed by the mother (with a coefficient of variation of 0.3). Simulations are performed as in the previous section, and performance colormaps can be computed for all possible concentrations of different ligand combinations (Fig. 6). Comparison of Fig. 6A-F with the corresponding Fig. 5A-F evidences that selectivity is more difficult to achieve in conditions of heritability (i.e., regions of optimal selectivity (marked in dark red) are reduced and occur at higher concentrations). Numerical values for the minimal concentration for the threshold of selectivity, as well as the reduction in total final concentration due to the combinatorial treatment in conditions of heritability, are shown in Table 3. Fig. 6 G-H, plots the values of AER and TER of the cells in the population before (Fig. 6G) and after (Fig. 6H) 60 hours of treatment for the minimal concentration of Ch 2 that meets the threshold in conditions of heritability. Comparison of both distributions with conditions of no heritability (Fig. 2J,M) evidences that heritability increases the variability in the expression levels of AER and TER in the population, resulting in a decrease in the performance of the drug combinations and a reduction in the region of optimal selectivity (Fig. 6A-F).

Conclusions and Discussion
Chimeric ligands with selective potential constitute one of the forefronts in modern pharmacology. The development of strategies to affect only malfunctioning cells inside a healthy tissue based on a sequential mechanisms of targeting is still in its early stages. Rational approaches based on modulating the strength of the interaction between ligand and target has shown that selectivity can be improved in a rational predictive manner [15,25,34]. Unfortunately, this results in a marked increase of the total concentration of drug that needs to be administered, which potentially increases the risk of toxicity and other undesired effects. Therefore, the problem of achieving selectivity at reduced drug concentrations is a main concern when developing selective drugs.
Our previous modeling approaches [25,34] allow us to predict the optimal value of the affinity and dissociation rates of both AE and TE for improved selectivity at the lowest drug concentration. Unfortunately, the affinity and dissociation rates in a given ligand-receptor interaction cannot be modulated gradually, since single mutations in the ligand change abruptly the binding and unbinding rates with the complementary receptor. In this sense, combination of two ligands can, in principle, result beneficial to improve the selective potential of the treatment since, for instance, highly potent ligands could affect cells expressing high TER concentrations, while more selective chimeras (i.e., with reduced potency in the AE subunit) could discriminate better between healthy and unhealthy levels of the target protein.
To our knowledge, our results constitute the first studies focused on the combination of selective drugs, by generalizing our previous results of single treatments with selective drugs [25] to study selective drug combinations in cell population models. Our studies show that the combination of selective drugs is synergistic in terms of their selective potential, i.e., the combination of selective drugs can reduce the total drug administered to achieve a given selective effect, compared to the same drugs acting alone. We also show that using a explicit model of two selective drugs is equivalent to a simplified model where the two drugs are assumed to interact additively. This alternative method allows us to develop performance maps where selectivity is computed for any given concentration of the combined drugs. We used a cell population based model to study how these types of treatments respond in a context of cell-variability and their robustness in condition where the amount of target proteins is inherited from mother to daughter cells. Interestingly, despite assuming additive interaction of the different chimeras when combined (i.e., they compete for the same molecular targets), when looking at the selective potential of the treatment, chimeras behave as synergistic.
Several main assumptions are taken into account while developing the model. First, we assume that production and degradation of each receptor is balanced in conditions of no ligand stimulation. We also simplified all potential downstream regulation in receptor expression after activation, focusing only on the regulation that takes place due to direct ligand stimulation. We also assume that the activity triggered by the ligand-receptor interaction is proportional to the amount of maximum active complexes formed. Other potential values such as the total value of active complexes at a given time also produce equivalent results, as discussed in [25] at the single cell level. Regarding the simulations of the population dynamics, we assumed that all cells proliferate at the same mean rate, independently of the amount of EGFR receptors. It is well-known that EGFR stimulation is correlated with the activation of proliferative signals [32], but experimental data monitoring differences in cell cycle length for Daudi versus Daudi-EGFR cells used to inform our model are not available [15]. In addition, the effect of heritability in the expression of receptors was assumed to simply depend on the amount of receptors expressed by the mother cells. Other potential scenarios to capture the effect of mutations in the regulation of the expression of receptors will be more realistic, but they will result in more free parameters and assumptions. In addition, we assume that the effect of heritability will be more relevant in longer experiments, i.e., when more generations of cells are allowed to develop. Unfortunately, experimental data are only available at the time point of t = 60h, corresponding to an average of 2.2 generations, insufficient to observe the selective pressure effect induced by the drug treatment. To mimic cell-to-cell variability in the population, we assumed gamma distributions for the amount of receptors expressed and the cell cycle length, based on several publications. Other types of distributions were also tested (gaussian, lognormal), with almost no difference in the results compared to the gamma distribution [35][36][37][38]. To quantitatively compare the different combinatorial treatments, a threshold is defined in terms of the potential selectivity of the treatment towards the different cell types expressing different concentrations of the target proteins (80% survival of the healthy cells while the number of unhealthy cells is maintained). Other potential threshold values defined also evidence the reported synergism when combining two selective ligands, but at different drug concentrations.
In conclusion, we have shown that combination of selective drugs can selectively affect a given cell population at reduced concentrations compared to single drug treatment. These types of theoretical studies focused on the rational design of selective drugs and treatments can complement experimental efforts, allowing researches to develop a more reliable and efficient approach to quantitative pharmacology.
Supporting Information S1 Text. Response surface plots for drug interaction in conditions of drug additivity, synergy or antagonism.
(TEX) S1 Fig. Workflow to obtain the effect of a drug combination on a cell population. (A) Direct simulation of two simultaneous treatments (see section Models: Eqs. 1-4 are solved directly for two simultaneous ligands (j = 1,2) at a constant concentration. The value of AER-AE complexes formed is then translated into a physiological effect using calibration with experimental dose response curves obtained from [15]. Two populations of cells are defined with values for AER and TER from gamma distributions for healthy and unhealthy cells. Eqs. 1-4 are solved numerically for each cell in the two populations, obtaining the dynamics of growth for healthy and unhealthy cell populations for a given constant concentration of L 1 and L 2 . (B) Calculation of the effect of combinatorial treatment assuming additive interaction between ligands: the maximum number of AER-AE complexes is calculated for each combination of AER and TER receptors concentrations by solving Eqs. 1-4 for a single ligand treatment (j = 1). The output of the model is translated to a calibration curve [25], obtaining the theoretical dose-response curves for each ligand. Physiological response curves are fitted to a four-parameter sigmoidal (Eq.9), and the physiological response for any concentration of two ligands is then calculated using the Loewe approximation for additive ligand interaction (Eq. 11). This response is then used to perform simulations for healthy and unhealthy cell populations, following the same procedure as in (A). Finally, the number of healthy and unhealthy cells after 60 hours of treatment is plotted in the corresponding isobologram for each ligand combination. The final performance colormap for each value of the combination of ligands is obtained by subtracting the normalized isobolograms for unhealthy minus healthy cells. Values above threshold of performance are highlighted in dark red.