The mathematical model

Background of the model

Doxo mediate its anti-tumor effect in a calcium-dependent manner. Doxo-induced remodeling of calcium signaling in TNBC has been observed previously [12, 13, 15, 20]. A schematic diagram (Figure 1) shows the calcium remodeling during the Doxo treatment in TNBC. Doxo treatment in TNBC increases the expression level of Ca²⁺ regulating proteins associated with calcium channels and pumps of the plasma membrane and ER namely, store-operated calcium entry (SOCE, ORAI1), TRPC1, IP₃R2, PMCA2, and SERCA1 [13]. The increased activity of ORAI1, which is the main calcium influx channel of the plasma membrane, is connected with therapeutic resistance, invasive phenotypes, and distant metastases in many cancers [24–27]. TRPC1 is a receptor gated calcium influx channel of the plasma membrane and is included in the proliferation of various cancer cell lines [28–30]. Another plasma membrane calcium channel is PMCA2 that vacillates the calcium from the cytosol into extracellular. Its increased activity has been observed in various studies during Doxo treatment [13, 31]. However, Doxo doesn’t significantly change activities of few other calcium channels namely, voltage-gated calcium channels (VGCCs) and Na⁺/Ca²⁺ exchanger (NCX). SERCA1 and IP₃R2 are ER influx and efflux calcium channels respectively. Their increased expressions have also been observed during Doxo treatment in studies done previously [13, 32].

Display full size

Figure 1. 

Schematic diagram: Doxo-induced modulation of calcium signaling pathways in TNBC cells [13]. The red arrow represents the increased mRNA level of calcium regulators during Doxo treatment

Analytical results

Result 1. (i) The solution of the system (1) with respect to the initial conditions (2) exists and is unique in the interval (0, ∞) and it remains positive for all future time t > 0 (see Supplementary A for proof); (ii) all solutions of system (1) with positive initial conditions are uniformly bounded within a region Γ, where Γ = {([Ca²⁺]_c, [Ca²⁺]_ER, [P_so]): 0 < [Ca²⁺]_c (t) <, $0<{\left[{\text{Ca}}^{2+}\right]}_{\text{ER}}\left(t\right)<\frac{\Lambda }{\gamma +\xi }$, $0<{\left[{\text{Ca}}^{2+}\right]}_{\text{ER}}\left(t\right)<\frac{p_2}{\delta }$, 0 < [ P_so] (t) < 1}, with $\Lambda =\alpha +\beta +\gamma c_{\mathit{es}}+P_2+\frac{p_1{p}_2}{2q_1\delta }$. Thus, the dynamics predicted are bounded within biological values (see Supplementary A for proof).

Result 2. The system (1) has interior equilibrium point $E^{*}\left({\left[{\text{Ca}}^{2+}\right]}_{\text{c}}^{*},{\left[{\text{Ca}}^{2+}\right]}_{\text{ER}}^{*},{\left[P_{\text{so}}\right]}^{*}\right)$ which can be obtained from the following expressions,

The interior equilibrium point exists if we obtain the positive value of ${\left[{\text{Ca}}^{2+}\right]}_{\text{c}}^{*}$ from the expression (3) and $p_1>\delta \left(q_1^2+{\left[{\text{Ca}}^{2+}\right]}_{\text{c}}^{*}{}^2\right)$ holds. It is locally asymptotically stable if a₂₁ > 0 and A₁A₂ > A₃ holds simultaneously (see Supplementary A for proof).

Choice of parameter set and validation of analytical results

A recent study done by Bong et al. [13] explored the role of Doxo through calcium signaling. They did a comparative analysis of cytosolic calcium dynamics (in MDA-MB-231 cell line) for three different conditions namely, control-treated cell line, 0.03 μM Doxo treated cell line, 1 μM Doxo treated cell line. They observed dynamics of cytosolic calcium using three criteria, which were time to reach the peak (peak time), peak amplitude, and AUC between control and Doxo treatments (0.03 μM, 1 μM). Altered calcium dynamics were observed in Doxo treated cells. The time to reach the peak and the peak amplitude of cytosolic calcium were observed between 5–20 s and 0.3–0.4 μM respectively in all three conditions. However, there was a significant difference in the AUC of control (AUC ≈ 22) and Doxo treatments (AUC ≈ 30). To study these altered calcium dynamics during Doxo treatment through our model (1), we simulated our model to replicate the above mentioned experimental results.

We explored different literatures and obtained few parameters while the rest of the parameters namely, β (SOCE calcium influx rate), γ ([Ca²⁺] influx rate through TRPC1 channel), δ(leaking rate), ξ([Ca²⁺] export through NCX), p₁ ([Ca²⁺] influx rate through IP₃R₂), q₁ (IP₃R₂ channel’s half saturation constant), p₂ (flux rate through SERCA1 channel) were estimated using experimental data points by Bong et al. [13]. Here we used least square method, which minimize the sum of the residuals to estimate the parameters. The base-line values of all the model parameters are given in Table 1. The model simulation results were also matched with the time to reach the peak, peak amplitude, peak width, and AUC as observed in the study [13], see Figure 2a. These cytosolic calcium dynamics are termed as control response dynamics (CRD). The dynamics of other state variables namely, [Ca²⁺]_ER and [P_SO] are also obtained and shown in Figure 2b. A histogram study was carried out to define the ranges of model readouts (peak amplitude, peak width, and AUC) for further computational exercises in our study. We varied the base set 10% up and down which represents the neighborhood of our baseline parameter set and for which we have observed no change in the dynamical output of the system [later we explored this in detail in the robustness exercise in section (Effect of different parameters on the system dynamics)]. The parameters were varied using the Latin hypercube sampling (LHS) technique that employs Monte Carlo random sampling technique. A total of 10,000 runs were done and for each run, we calculated the peak amplitude, peak width, and AUC. A histogram plot in Figure 2c shows the range of the above three readouts establishing the range of CRD. These ranges will be used in future simulations for capturing the CRD.

Display full size

Figure 2. 

Preliminary analysis: (a) the best fit curve obtained by the least square method for estimation of parameters using experimental data [13]; (b) the dynamics of all state variables of the model for the parameter set given in Table 1; (c) the histogram plots of peak amplitude, peak width, and AUC were obtained on the parametric variation of ± 10% from the baseline; (d) the 2D parameter space between β and p1. Here, green color represents stable (focus) CRD region, red color represents stable (node) CRD region with negative real eigenvalues and blue color represents non-CRD stable region. The black star mark in the figure represents the baseline value of both parameters

We also observed that for this baseline parameter set (Table 1) the existence conditions for the interior equilibrium point E* was satisfied, and the interior equilibrium point E* = (0.0903, 11.058, 0.9976) is locally asymptotically stable as the eigenvalues of the characteristic equation is (−36.5149, −0.0324 – 0.0033i, −0.0324 + 0.0033i). This follows that CRD shows stable dynamics through damped oscillations. SOCE is an important plasma membrane channel of calcium entry in the cytosol that opens in response to ER calcium depletion. Entry through this channel is mediated by STIM1 proteins that disseminate from the ER membrane and binds to ORAI1 proteins and activates SOCE. IP₃R2 channel plays an important role in the exit of calcium from ER into cytosol. The involvement of IP₃R2 channel in mediating SOCE has been established earlier thereby influencing cytosolic calcium dynamics [41, 42]. To observe this complex association between β channel’s influx rate) and p₁ ([Ca²⁺] influx rate through IP₃R2) and its role in establishing CRD. We varied these two parameters 2-fold up down from their base values and generated the 2D parameter space shown in Figure 2d. The green region represents that the system achieves CRD and stability through damped oscillations. Whereas, the system is a stable node in nature and attains CRD in the red region. We observed that on increasing β and setting p₁ to its baseline value, the model dynamics changes from stable focus to stable node dynamics.

Dose-dependent cytosolic calcium profiles during Doxo treatment

Doxo works in a concentration-dependent manner [13]. An altered calcium dynamics was observed depending on the concentration of the Doxo given in the MDA-MB-231 TNBC cell lines. They observed a significant difference in calcium dynamics of control and Doxo treated cells. Among the readouts defined in the earlier section, the most significant indicator differentiating the two cases was found to be AUC [13]. Higher AUC value in the case of Doxo represented a more sustained cytosolic calcium profile. To obtain Doxo induced calcium remodeling observed in the study done by Bong et al. [13], we varied the parameters related to Doxo associated proteins. Doxo remodels the calcium dynamics by increasing the activity/expression levels of a few calcium channels and pumps, hence increasing the calcium flux flow through these pumps and channels. To incorporate Doxo associated remodeling in our model, we also employed a similar technique to simulate such effects through our model. The parameters associated with Doxo are β (SOCE calcium influx rate), γ ([Ca²⁺] influx rate through TRPC1 channel), p1 ([Ca²⁺] influx rate through IP₃R2), p2 (flux rate through SERCA1 channel), p₃ (flux rate through PMCA2 channel), and they were simultaneously and step-wise increased by 25% (Doxo 1), 50% (Doxo 2), and 75% (Doxo 3). The calcium profile generated through elevated parameter ranges was given in Figure 3a. We found that the 75% increase in the parameters related to Doxo associated proteins reproduced the cytosolic calcium profile observed experimentally in Doxo treated cells by Bong et al. [13]. We also observed that AUC increases with every 25% step increase in the parameters related to Doxo associated proteins, see the bar graph [Figure 3b]. The AUC values for all four cases were in 90% confidence interval (CI) of the mean value (see Table 2). Moreover, the 90% CI in all four cases were non-overlapping showing that the observed trend of increase in AUC is also statistically significant (for details see Supplementary B).

Display full size

Figure 3. 

Dose dependent analysis: (a) dynamics of cytosolic calcium for four different parameter sets. The blue curve in the figure shows the cytosolic calcium dynamics obtained from baseline parameters while red, yellow, purple curves show the cytosolic calcium dynamics after increasing the parameters related to Doxo associated proteins from base values (as given in Table 1) by 25%, 50%, 75% respectively and peak values of cytosolic calcium (in μM) are also shown; (b) the bar graph shows the AUC values obtained in the four cases. Here Doxo 1, Doxo 2, Doxo 3 denotes the categories when the baseline parameter set is increased by 25%, 50%, 75% respectively

Effect of different parameters on the system dynamics

In this section, we explored the effect of parameter variations on the system’s dynamics. We first varied single parameters through robustness analysis and observed their effect, and then we performed a global sensitivity analysis (GSA) to capture the effect of multi-parameter variation.

Robustness analysis

The robustness analysis was done to assess the model’s dynamics to parameter perturbations and to establish the relevance of the obtained parameter set. To investigate whether the parametric variations disrupt the system dynamics or not, each parameter was increased and decreased one at a time by the step size of 1% of its baseline value up till 100% (see Figure 4). The bar length indicates the robustness of the parameter in conserving the model dynamics obtained from our baseline parameter set. A lower bar length specifies that a proportionately small perturbation hampered the model dynamics, while a bigger bar length specifies the robustness of the parameter in retaining the dynamical nature of the model. Our model was robustly producing the CRD against the parametric perturbations except in the case of parameters ξ, p1, q₃ (see Table 3). Here, ξ represents the export of calcium into extracellular through NCX channel, p₁ represents the calcium flux rate of IP₃R2 channel and q₃ represents half-saturation constant of PMCA2 channel.

Display full size

Figure 4. 

Robustness analysis: the figure shows the dependence of the model’s dynamics on parameter variations. Parameters were varied 100% up down from the baseline values. The length of bars represents the range of parameters up to which the model dynamics hold despite parameter variation. The green color bar represents the robustness result achieved by increasing the parameter value and the yellow color for the opposite case

Table 3.

Robustness analysis table: effect of individual parameter variation on dynamics of the whole system. Here CRD and non-CRD indicate control response dynamics and non-control response dynamics respectively. By varying each parameter value individually, we observed and recorded the ranges for CRD model dynamics and non-CRD model dynamics

Parameter	Range	Dynamics
α	0 < α < 0.08	stable CRD
β	0 < β < 10.40	stable CRD
γ	0 < γ < 0.0002	stable CRD
δ	0 < δ < 0.0028	stable CRD
ξ	0 < ξ < 1.33	stable non-CRD
	1.33 < ξ < 9.8	stable CRD
	9.8 < ξ < 14	stable non-CRD
p1	0 < p1 < 0.009	stable non-CRD
	0.009 < p1 < 0.031	stable CRD
	0.031 < p1 < 0.035	stable non-CRD
q1	0 < q1 < 0.1036	stable CRD
p2	0 < p2 < 7.2336	stable CRD
q2	0 < q2 < 0.2	stable CRD
p3	0 < p3 < 14.68	stable CRD
q3	0 < q3 < 0.6	stable non-CRD
	0.6 < q3 < 3	stable CRD
ces	0 < ces < 2000.3624	stable CRD
ks	0 < ks < 100	stable CRD
τs	0 < τs < 60.3504	stable CRD

Display full size

Global Sensitivity Analysis

The GSA identifies the efficiency of model parameters and delivers crucial details about the system performance. It encapsulates the parameters whose perturbation significantly affects/correlates the model output dynamics. Here we chose partial rank correlation coefficient (PRCC) method, which is based on sampling, for GSA. In our study, we used the LHS method as it compactly differentiates the input parameters [43]. PRCC measures the association between model inputs and outputs using correlation through the LHS method [44–46]. Each parameter was sampled over a set of values obtained by varying parameters ± 10 fold of their baseline value (Table 1). Using LHS, 10,000 sets of parameters were generated and PRCC was evaluated varying between −1 and 1. A cut off ± 0.3 was set to the sensitivity index acquired above for refining sensitive parameters [47, 48].

We performed sensitivity analysis for all state variables of our model (see Figure 5a). We observed that the SOCE influx rate (β), maximum Ca²⁺ influx through TRPC1 channel (γ), Ca²⁺ export through NCX (ξ) and calcium concentration in extracellular space (ces) were significantly correlated with cytosolic calcium (see Table 4). Whereas, maximum Ca²⁺ influx through TRPC1 channel (γ), Ca²⁺ export through NCX (ξ) and calcium flux rate (p1) and half saturation constant (q₁) of IP₃R2 channel, calcium flux rate (p2) and half saturation constant (q₂) of SERCA1 channel were significantly correlated with [Ca²⁺]_ER (see Table 4).

Display full size

Figure 5. 

(a) GSA: the PRCC index plotted at t = 200 s shows the sensitivity of each model parameter against the state variables of our model. The higher the bar length, the more sensitive is the parameter to that variable. The parameter correlated negatively as well as positively with the model state variables; (b) temporal sensitivity analysis: the figure shows the temporal sensitivity analysis of parameters associated with calcium channels and pumps whose activities or expressions are significantly increased during Doxo treatment [13]

Table 4.

Global uncertainty analysis results: the table shows the sensitivity indexes of model parameters for all three state variables

Parameters	[Ca2+]c	[Ca2+]ER	[PSO]
α	0.0839	0.0341	−0.0048
β	0.5562*	0.2850	−0.0420
γ	0.6332*	0.3585*	−0.1129
δ	0.0230	−0.0678	0.0409
ξ	−0.8523*	−0.6793*	0.2704
p1	0.0986	−0.5638*	0.2487
q1	−0.0904	0.5954*	−0.1636
p2	−0.0793	0.5942*	−0.2100
q2	0.1606	−0.7296*	0.2751
p3	−0.0458	−0.0206	0.0217
q3	0.1455	0.0842	−0.0390
ces	0.6351*	0.3520	−0.0873
ks	0.0957	0.0277	0.3719*
τs	−0.2459	−0.1178	−0.8330*

Display full size

Each parameter was sampled for 10 fold up and down from the baseline value and obtained sensitivity indexes using PRCC algorithm. A cut-off ± 0.3 was applied over sensitivity indexes and were marked sensitive [47, 48].

*

: represents the value beyond the cut-off ± 0.3

It is known that Doxo significantly increases the mRNA levels of calcium regulators associated with ORAI1, TRPC1, IP₃R2, SERCA1, and PMCA2 [13]. Temporal sensitivity analysis was performed to study the temporal effect of the parameters (β, γ, p₁, p2, p₃) associated with these calcium channels and pumps in our model (see Figure 5b). We found that the parameters had a time-dependent sensitivity towards cytosolic calcium. We observed that two parameters β and p₁ exchanged their sensitivity with time. Initially, p1 was sensitive but its sensitivity reduced with time, whereas the sensitivity of β, γ increased with time. The parameter p₂ lost its sensitivity very fast after the first 10 s, while the parameter p₃ didn’t change its sensitivity and remained insensitive throughout the time course. These results were suggesting the varying role of different calcium pumps in remodeling cytosolic calcium profile.

Parameter recalibration

The remodeling in calcium dynamics caused by Doxo treatment in TNBC significantly increases the activity of calcium pumps and channels. This resulted in sustained cytosolic calcium levels thereby increasing the area under the curve of cytosolic calcium (see Figure 3b). This sustained calcium level in cytosol has been known to cause oncogenesis in previous reports [14, 15]. Moreover, remodeling of Ca²⁺ channels and pumps regulating proteins are connected with drug resistance in cancer [19, 49, 50]. Wang et al. [49] showed that increased activities of TRPC5 due to which multi-drug resistance pump, P-glycoprotein up-regulates in MCF-7 breast cancer cells resulted in Doxo resistance. The study showed that suppression of TRPC5 calcium channel reversed Doxo resistance in tumors formed in the MCF-7 cell line. Therefore, to mitigate the calcium-associated block in Doxo treatment, we performed parametric recalibration aiming to find parameters whose perturbations can restore CRD. We took the parameter sets which produced cytosolic calcium profiles obtained in section (Dose dependent cytosolic calcium profiles during Doxo treatment), where an increased area under the curve was observed when parameters related to Doxo associated proteins were perturbed by different step size. The range of readouts (peak amplitude, peak width, and AUC) from histogram plot (see Figure 2c) were 0.28–0.34 μM, 24–28 s and 23–27, respectively. These readout ranges are termed as CRD ranges and were further used as matching criteria for the parametric recalibration, i.e. on varying parameters, if the profile of cytosolic calcium so obtained had the above-mentioned attribute values, then it was declared that the CRD has been achieved.

We performed a single parameter variation where each parameter was perturbed 100% up-down from their respective values and we filtered out parameters that retained the Doxo profile to the control profile (CRD). This exercise was repeated for different Doxo treatments where associated parameters were increased by 25% (Doxo 1), 50% (Doxo 2), 75% (Doxo 3). Single parametric variation was able to restore the CRD only in the case of 25% and 50% increase. The results of single parameter variation for the first two cases of Doxo treatments are given in Figure 6(a–b). For parametric recalibration in the 75% case, we used two-parameter variation technique. For this rather going by brute force approach, we chose parameters related to Doxo associated proteins of our model that were, β (SOCE influx rate), γ (calcium influx rate of TRPC1), p1 (IP₃R2 calcium influx rate), p₂ (calcium influx rate of SERCA1), p₃ (PMCA2 calcium influx rate). We picked one parameter from the above and reduced its value by 20%. After that, all other parameters were varied 100% up and down from their respective values (i.e. in the Doxo 3 parameter set) as done in the single parameter variation. Out of the five parameters related to Doxo associated proteins, the reduction of β (SOCE calcium influx rate) by 20% turned out to be the only parameter that can restore CRD in combination with other model parameters. Except for β (SOCE calcium influx rate), reduction in other parameters related to Doxo associated proteins was not sufficient to restore CRD in combination with any other model parameter. The results of double parameter variation in 75% increased case are given in Figure 6(c).

Display full size

Figure 6. 

Parametric recalibration results: (a) the figure plot shows the single parameter variation results for Doxo 1; (b) the figure plot shows the single parameter variation results for Doxo 2; (c) the figure plot shows the double parameter variation results for Doxo 3 after reducing β (SOCE calcium influx rate) by 20%. The green bar shows the range for non-CRD and the yellow bar shows the range after CRD is obtained

We filtered out those restoration approaches which robustly produced CRD once obtained in our parametric variation range with 90% CI of readout values (for details see Supplementary B) in the range of CRD and non-overlapping CI from their respective initial parametric conditions (Doxo 1, Doxo 2, Doxo 3). The consolidated results are shown in Table 5.

In silico screening of calcium related drugs

Breast cancer cells display innate and acquired resistance to numerous anticancer drugs [51], which is the main obstacle in the effective treatment of breast cancer [52]. Therefore, enhanced medication protocols and alternative chemotherapeutic approaches are required. Studies done previously, observed the significant synergistic effects of simvastatin and Doxo together in breast cancer cells and suggested that combining simvastatin and Doxo together can be an effective treatment for breast cancer [20]. This calcium-mediated overcoming Doxo resistance encouraged us to screen known calcium-related drugs that can be administered with Doxo.

To execute this in silico experiment, we employed the CRD restoration approach (i.e. the three readout values are in the range of CRD) used in the recalibration exercise done in the previous section (Parameter recalibration). We downloaded a list of drugs with their corresponding drug targets from Drug Bank database [53] that includes approximately 13,000 drugs and their corresponding targets. We filtered out 47 drugs (inhibitors) whose targets involve Doxo associated calcium channels and pumps present in our model namely, IP3R2, SERCA1, VGCC, ORAI1, TRPC1, and NCX. We also did the literature survey and found inhibitors of calcium channels and pumps present in our study [54–56]. A complete list of drugs, drug targets, and associated parameters of our model used for further analysis is given in the Supplementary Table. We defined two cases Doxo 3 and Doxo 3 + drug in which the associated parameter with the drug is reduced by 50%. The model readouts for 47 drugs in combination with Doxo were evaluated and then examined whether CRD is restored back or not. Out of 47 drugs, Pyr6 (pyrazole compound), linoleic acid, SKF96365 showed complete restoration of CRD, i.e. corresponding to these drugs, the readout values were in the range of CRD. Dantrolene and caffeine turned out to be the only drugs with peak width and AUC values in the range of CRD. The in silico screening results are shown in color-coded representation (see Figure 7). The yellow color represents that the particular readout value is in the range of CRD for the corresponding drug whereas the blue color represents that the particular readout value is not in the range of CRD.

Display full size

Figure 7. 

In silico drug screening: 47 calcium-related drugs were screened on the basis of the CRD restoration approach built-in the previous section. Three readouts namely, peak amplitude, peak width, and AUC were evaluated for each drug. The yellow color for a particular readout shows that it is in the range of CRD for the Doxo 3 + corresponding drug combination therapy and the blue color shows that it is not in the range of CRD

The 90% CI of readout values (see Supplementary B) was evaluated for five parametric conditions, i.e. the combination of Doxo 3 and drugs (Pyr6, linoleic acid, SKF96365, dantrolene, and caffeine). For Pyr6, linoleic acid, and SKF96365, the CI of readout values were in the range of CRD whereas, for dantrolene and caffeine, 90% CI of peak width and AUC were in the range of CRD. Moreover, in presence of the drug the CI of readouts was not overlapping with CI of readouts in the case of Doxo 3 implying a significant difference between Doxo 3 and Doxo 3 + drug (see Table 6).

Understanding effect of Doxo on calcium dynamics through mathematical model

In the present study, we proposed a mathematical model to capture the alteration in calcium dynamics due to Doxo treatments and predict possible restoration strategies. Our model was able to replicate the experimentally observed cytosolic calcium dynamics both in control and Doxo treatments. Initially, the calcium peak was observed around 2 s and later the system showed stability through damped oscillations in the control profile that was ascertained by our analytical results. The initial peak amplitude is associated with enhanced SOCE activity and store depletion [15, 57]. The results from robustness analysis showed that our system was able to reproduce CRD against parametric perturbations except for ξ, p1, q₃. The NCX exchanger together with calcium ATPase pumps plays an important role in restoring resting cytosolic calcium [15]. We defined the robust ranges for each model parameter in reproducing CRD. Moreover, through ± 10% parametric perturbation, we revealed the range of CRD in the robust parametric neighborhood.

Doxo works in a dose-dependent manner showing anti-proliferative effect at lower concentration and inducing cell apoptosis at higher doses [58, 59]. However, higher dosage of Doxo results in calcium remodeling in breast cancer cells [13, 19–21]. The Doxo associated calcium pumps and channels have taken in our model were significantly perturbed in higher dosage as compared to the lower dosage of Doxo [13]. Prolonged treatment of Doxo at higher dosage also results in overexpression of STIM/ORAI channel that has known to show anti-apoptotic effects [60]. A dose-dependent study of Doxo done by Bong et al. [13] revealed different modes of cell death based on the cytosolic calcium dynamics obtained at different Doxo doses. At higher dosage, Doxo shows its apoptotic effect in a calcium dependent manner hence this dose-dependent remodeling of calcium can affect the overall cellular fate. However, short-term and long-term effect of Doxo remains to be completely elucidated. We tried to mimic the experimentally obtained dose-dependent calcium remodeling by Doxo [13]. We evaluated AUC for every dose-dependent curve in our study. We observed that on increasing Doxo dose, i.e. on increasing the parameters related to Doxo associated proteins, AUC increased. This implies a positive correlation between Doxo doses and AUC. The dose-dependent curve obtained on increasing parameters related to Doxo associated proteins by 75%, completely replicated the curve obtained experimentally by Bong et al. [13]. The increased/prolonged cytosolic calcium during treatment is the biological relevance of increased AUC on increasing parameters related to Doxo associated proteins by 75%. The prolonged elevation of cytosolic calcium can influence the ability of cells to monitor progression through the cell cycle and lead to unchecked proliferation and tumorigenesis [61]. This motivates us to explore more on understanding the influence of parameters related to Doxo associated proteins on cytosolic calcium dynamics.

Role of Doxo associated calcium channels and pumps

The activity of various calcium channels and pumps increases during Doxo treatment [13]. Many of these channels and pumps are associated with drug resistance [17–19]. The results of sensitivity analysis showed a high sensitivity of SOCE influx rate (β) and TRPC1 influx rate (γ) on cytosolic calcium. Moreover, these channels influx rates are among parameters related to Doxo associated proteins. The channel SOCE is an important plasma membrane calcium channel that activates in response to ER calcium depletion. SOCE is involved in various oncogenic processes [24], e.g., cancer cells invasion in glioblastoma and increase in distant metastasis in breast cancer [25]. Abdoul-Azize et al. [20] reported an increase in cytosolic calcium levels in MDA-MB-231 breast cancer cells after Doxo treatment due to increased activity of STIM1 and TRPC1. This underlines the crucial role of these channels in establishing a transient nature of cytosolic calcium during Doxo treatment. The extracellular calcium concentration also turned out to be sensitive to cytosolic calcium. This may be because, during Doxo treatment, MDA-MB-231 cells showed an enhanced sensitivity towards extracellular calcium perturbations [13]. The sensitivity analysis results portrait the calcium remodeling during Doxo treatment. Further, the temporal sensitivity analysis revealed the temporal landscape of parameters related to Doxo associated proteins for cytosolic calcium. The heterogeneity in the responses against the parameters related to Doxo associated protein perturbations at early and later time points elucidated the specific role of each calcium channel/pump on cytosolic calcium dynamics. The rise in the sensitivity of SOCE calcium influx rate (β) and TRPC1 calcium influx rate (γ) at later time points could be because they are calcium influx rates thereby, increasing cytosolic calcium levels during Doxo treatment. Interestingly, the sensitivity of ER related calcium pumps with flux rates of p₁ and p₂ decreased with time, and these two parameters were significantly over-expressed in Doxo treated cells [13]. The magnitude of calcium release from the ER calcium store into cytosol depends on the difference between calcium concentration of ER and cytosol [62]. Due to high influx rates from store-operated channels, the cytosolic calcium levels are elevated in Doxo treated cells. As a result, at later time points, there is a decrease in calcium concentration gradient between ER and cytosol, i.e. ER calcium levels decrease with time in correspondence to cytosolic calcium. Hence, the ER flux rates p₁ and p₂ become less sensitive at later time points. Each of these calcium channels and pumps modulated by Doxo holds a specific function in modulating calcium dynamics during Doxo treatment. Thus, targeting these dysregulated Ca²⁺ channels and pumps might provide efficient chemotherapy in TNBC.

Restoration of CRD and indispensable role of SOCE (ORAI1) calcium influx rate

A major goal of the present study was to unveil possible mechanisms that can overcome the calcium-mediated Doxo drug resistance. This might lead to the identification of novel therapies. We hypothesized that restoration of calcium CRD can reduce the Doxo drug resistance. Our dose-dependent restoration approach showed diverse treatment strategies. The treatment strategies which we found relevant at lower Doxo dose were not sufficient to restore CRD at higher Doxo doses. The single parameter variation was adequate to restore back CRD for moderate cases (where parameters related to Doxo associated proteins are increased by 25% and 50%). The common restoration technique in these two cases involves a decrease in the flux through calcium channel TRPC1 (γ). The role of TRPC channels in calcium-related drug resistance has been reported previously where a block in these channels increased the cell sensitivity towards Doxo treatment [19]. At a higher Doxo dose, i.e. 75%, we found the indispensable role of SOCE flux rate parameter (β). It was only after the decrement in the value of β by 20%, parametric combinations were able to restore the CRD. Hence, β came out as the primary target in the path of overcoming calcium-mediated Doxo drug resistance. The results from the sensitivity analysis also emphasized the high calcium sensitivity against β perturbation and its role in cancer proliferation has been explored in studies earlier [57]. Hence, the indispensability of SOCE calcium influx rate (β) is paramount in the Doxo drug resistance.

Identification of compounds as potential anti-cancer drugs in combination with Doxo

Doxo is a broad-spectrum chemotherapeutic drug used in the clinic [63]. Targeting the altered activities of Ca²⁺ channels and pumps may present promising chemotherapy for cancer treatment. The prominent research towards the assessment of various Ca²⁺ channel blockers, regulators, and inhibitors as anti-cancer drugs is emerging nowadays [14]. The in silico drug screening methodology employed in our study identified the drugs which together with Doxo can overcome calcium-dependent Doxo drug resistance. Calcium compounds Pyr6 (pyrazole compound), linoleic acid and SKF96365 were the only successful candidate which restored CRD. Linoleic acid is given as nutritional supplementation and for treating dietary deficiency or irregularity [64]. In a recent study, pyrazole derivatives have been considered as selective blockers of TRPC channels [65]. Pyr6 is a pyrazole compound that inhibits ORAI1 mediated calcium entry into the cytosol (represented by β in our model). Linoleic acid selectively inhibits the activity of the TRPC1 channel (represented by γ in our model) but studies have shown that it induces some STIM1-STIM1 association while inhibiting STIM1-ORAI1 binding [53, 66]. The results obtained in the drug screening exercise were in concordance with the results obtained from the recalibration exercise. The putative drug targets of the successful candidates were similar to the suggested mechanism for restoration. Linoleic acid has been screened in case of breast cancer. It has shown anti-tumor properties in breast cancer cells [67]. SKF96365 has already been known to have anti-tumor properties in a calcium-dependent manner. On treating MDA-MB-231 human breast tumor cells with SKF96365 (which blocks the store-operated [Ca²⁺] influx), the proliferation of invasive tumor cells was reduced [25]. Moreover, SKF96365 has been known to attenuate drug resistance in combination with Doxo in hepatocellular carcinoma cells [68]. This underlines the importance of this combination therapy as a putative strategy in overcoming drug-resistance.

The present study focused on the mechanistic viewpoint of reduced sensitivity of TNBC towards Doxo. The investigation highlights the importance of integrating the calcium signaling of various calcium regulating compounds for their effective anti-tumor effects deliverance along with chemotherapeutic agents. Here, we proposed few potential strategies which can be used in combination with Doxo. We further used a drug dataset to identify three drugs from existing literature that could be used individually along with Doxo as a combination therapy. Before ending this article, we would like to mention that our work lacks experimental validation, but the proposed approach delineates the complex nature of Doxo drug-resistance and can help in developing effective therapeutic strategies for TNBC.