The Quest to Quantify Selective and Synergistic Effects of Plasma for Cancer Treatment: Insights from Mathematical Modeling.

Cold atmospheric plasma (CAP) and plasma-treated liquids (PTLs) have recently become a promising option for cancer treatment, but the underlying mechanisms of the anti-cancer effect are still to a large extent unknown. Although hydrogen peroxide (H2O2) has been recognized as the major anti-cancer agent of PTL and may enable selectivity in a certain concentration regime, the co-existence of nitrite can create a synergistic effect. We develop a mathematical model to describe the key species and features of the cellular response toward PTL. From the numerical solutions, we define a number of dependent variables, which represent feasible measures to quantify cell susceptibility in terms of the H2O2 membrane diffusion rate constant and the intracellular catalase concentration. For each of these dependent variables, we investigate the regimes of selective versus non-selective, and of synergistic versus non-synergistic effect to evaluate their potential role as a measure of cell susceptibility. Our results suggest that the maximal intracellular H2O2 concentration, which in the selective regime is almost four times greater for the most susceptible cells compared to the most resistant cells, could be used to quantify the cell susceptibility toward exogenous H2O2. We believe our theoretical approach brings novelty to the field of plasma oncology, and more broadly, to the field of redox biology, by proposing new ways to quantify the selective and synergistic anti-cancer effect of PTL in terms of inherent cell features.

1. Introduction

In the last decade, the use of cold atmospheric plasma (CAP), which is an ionized gas near room temperature, has become a novel method to treat cancer. Both direct application of CAP (e.g., by the clinically approved kINPenMED® plasma jet) and indirect treatment by application of plasma-treated liquids (PTLs) have been shown to provide a significant anti-cancer effect [1]. Van Boxem et al. [2] showed that PTLs have an anti-cancer effect for a number of different CAP and liquid conditions, and Lin et al. [3] found that CAP can induce immunogenic cancer cell death. This mode of cell death induced by CAP was later attributed to the CAP generated short-lived reactive species [4]. Moreover, CAP and PTLs have been reported to cause a selective anti-cancer effect [5], although selectivity depends on the cell type, the type of cancer, and the culturing medium [6]. Bekeschus et al. [7] demonstrated, using an in ovo model, that CAP is a safe cancer treatment modality with respect to possible metastasis formation. A number of promising results of the clinical application of CAP for cancer treatment have also been published (see e.g., [8,9]).

It is widely believed that the processes leading to cancer cell death are initiated by reactive oxygen and/or nitrogen species (RONS), in particular hydrogen peroxide (), but the knowledge about the specific mechanisms underlying cell death induced by CAP and PTL is still very limited. The lack of understanding of the combined effect of RONS contained in CAP and PTL in terms of the cellular response to exposure is problematic in the development of CAP/PTL treatment as a standardized cancer therapy for clinical use. Ultimately, it should be possible to predict and quantify the susceptibility to CAP/PTL of a particular cell line in terms of features specifically associated with those cells.

So far, the vast majority of the literature in plasma oncology are experimental studies. As a complement, another approach to increase the understanding of complex biological systems such as the interaction between cells and PTL, is to develop a mathematical model that includes all the known information (of major importance in the given context) about the system and use it to investigate the system’s response to various conditions. Especially, the system’s response to a perturbation of the “normal” conditions can be analyzed. Furthermore, the development of the mathematical model itself can be seen as a way to summarize the current state of knowledge on the matter in a compact manner; it can be seen as the current “working hypothesis” of the mechanisms and processes governing the system dynamics.

Mathematical modeling has indeed proven to be a useful approach to increase our knowledge about the mechanisms of the cell’s antioxidant defense and redox signaling. Some examples are the range of diffusion of in the cytosol [10,11] and the cellular decomposition of exogenous [12,13,14]. In the context of plasma oncology, two catalase-dependent apoptotic pathways associated with cancer cells, which possibly could be reactivated by CAP and thus explain the anti-cancer effect of CAP, have been investigated by mathematical modeling [15]. It was found that these pathways are unlikely to account for the anti-cancer effect of CAP and thus the underlying cause has to be studied further.

In the present study, we develop a mathematical model that includes the species and mechanisms of major importance in the context of a cell system exposed to PTL. The ultimate aim is to find a measure in terms of key features and characteristics of cells, which is able to quantify a particular cell system’s susceptibility toward PTL and thus explain differences in response between normal cells and cancer cells. To the best of our knowledge, this is a completely novel approach in the field of plasma oncology, and we believe that our study will provide a new perspective and new insights as a complement to the experimental studies. An extensive summary of the background, leading to the more detailed research question, is provided in Section 2.

2. Experimentally Observed Cytotoxic Effects of CAP and PTL and Possible Features Determining Cancer Cell Susceptibility

An immediate effect of CAP treatment of cancer cells is an increase in intracellular RONS [16,17,18,19,20]. The significance of this RONS accumulation has been verified by the observation that the treatment does not succeed if the cancer cells have been pre-treated with intracellular RONS scavengers [17,21,22]. The origin of the increase in intracellular RONS after CAP treatment is still under investigation, but a hypothesis consistent with experimental observations is that it is caused by a diffusion of extracellular CAP-originated RONS across the cell membrane [17,18,23,24].

It has been demonstrated that the anti-cancer effect of CAP can also be induced by the species in PTL. In PTL, which mainly consists of , and [25,26,27], has been shown to be of major importance [2,25,26,27,28,29,30,31]. It has been demonstrated that the consumption rate, which is cell specific, of cancer cells after PTL treatment, is a key factor determining the specific susceptibility of cancer cell lines to PTL. More explicitly, it has been reported that the higher the consumption rate of cancer cells, the lower the susceptibility toward CAP/PTL [32]. The susceptibility of cancer cells toward exogenous has also been shown in [33,34,35,36]. However, it has been found that alone cannot account for the total anti-cancer effect observed for PTL [27]. In this context, there are some reports of a synergistic effect of and in PTL [25,26]. Thus, the cytotoxic effect of seems to be enhanced in the presence of . The study in [26] found a selective, synergistic anti-cancer effect for in the μM-range and in the mM-range, whereas in [25], a non-selective, synergistic anti-cancer effect was reported when and were both in the mM-range. Since and in PTL may react to form [37], which is known to be highly toxic to cells, it has been speculated whether is the species causing the synergistic effect of and . Its formation could thus potentially increase the cytotoxicity of PTL compared to an equal concentration of only.

To summarize, some key points of the observed cytotoxic effects of CAP or PTL are:

An intracellular increase of RONS, which is likely to be caused by diffusion of CAP-originated constituents through the cell membrane, is crucial for cell cytotoxicity.

The key species in the anti-cancer effect of PTL is (note that this may be different for direct CAP treatment, where short-lived RONS also play a crucial role [4,38]), and the corresponding cytotoxicity is inversely proportional to the extracellular consumption rate of .

The effect of extracellular is enhanced in the presence of , which can be a clue to understanding why PTL enables a more efficient treatment than a mock solution of only.

Thus, from the information presented in literature, we can conclude that the cellular response to an addition of extracellular , with and without a simultaneous addition of , is crucial to understanding the anti-cancer effect of PTL. In Section 2.1, we introduce the key parameters to predict the response toward extracellular of cells. We also relate this to general differences between normal cells and cancer cells. Based on this information and knowledge, we introduce our approach, and formulate our research question and aim in detail in Section 2.2.

2.1. Differences in Cellular Response to Exogenous Hydrogen Peroxide

In particular, two factors determine whether a certain cell line is susceptible to the exposure of exogenous :

The plasma-membrane diffusion rate constant and

The intracellular expression of catalase.

Several cancer cell lines have shown a common phenotype of decreased catalase expression and increased aquaporin expression (that facilitates the transport of through the cell membrane [39,40,41,42] and thus determines the membrane diffusion rate) compared to normal cells. Hence, cancer cells in general can be assumed to be more susceptible to exogenous .

2.1.1. Membrane Diffusion Rate of Hydrogen Peroxide in Normal versus Cancer Cells

Aquaporins are proteins that form pores in the cell membrane. Primarily, they facilitate the transport of water between cells, but they also enable the trans-membrane diffusion of (due to the chemical similarities between both molecules). Thus, the aquaporin expression in the cell membrane relates to the membrane diffusion rate of . Many aquaporins have been found to be overexpressed in tumors of different origins, especially in aggressive tumors [43]. Since different cancer cell lines express aquaporins to various extents [43,44], the different responses of -exposure by different cancer cell lines can at least partly be explained by the non-identical levels of aquaporin expression. In [45], it was found that aquaporin 3 accounted for nearly 80% of the membrane diffusion of in a human pancreatic cancer cell line. For cells with a decreased aquaporin 3 expression, the rate of -uptake from the extracellular compartment was significantly decreased. It has furthermore been shown that for glioblastoma tumor cells, the anti-cancer effect of PTL as well as the increase of the concentration of intracellular RONS was significantly inhibited when aquaporin 8 was inhibited [46].

2.1.2. Catalase Activity in Normal versus Cancer Cells

Catalase is one of the main enzymes of the antioxidant defense system of cells of almost all aerobic organisms. The biological role of catalase is to regulate intracellular steady-state concentrations of , and experimental investigations and kinetic models using in vitro data have demonstrated that catalase is the major enzyme involved in the antioxidant defense against high concentrations of [12,47,48,49]. In particular, catalase has been shown to be responsible for the clearance of exogenous  in vitro and in vivo [12,50,51,52].

Although catalase levels vary widely across cell lines, the total concentration of catalase (extracellular and intracellular) is frequently reported to be lower in cancer cells than in normal cells [36,53,54,55,56,57,58,59,60]. In [61], it was found that the catalase activity in various cancer cells is up to an order of magnitude lower compared to normal cells, and in [62], it was shown that normal cells had a better capacity to remove extracellular than cancer cells; the rate constants for removal of extracellular were on average two times higher in normal cells than in cancer cells. Furthermore, it was reported in [62] that the rate constants for removal by different cell lines correlated with the number of active catalase monomers per cell.

However, while in general, the levels of catalase are low in cancer cells, catalase activity appears to vary greatly across different cancer cell lines [63]. In [34], it was found that three cancer cell lines (glioblastoma) that were extremely susceptible to (generated by ascorbic acid) had reduced activity of intracellular catalase. Ascorbic acid-resistant cancer cell lines, on the other hand, exhibited significantly higher levels of catalase, but catalase knockdown sensitized these cell lines to extracellular .

An additional aspect of catalase that may be of interest in the context of cytotoxicity of CAP and PTL, is that it has been shown to decompose [64]. Thus, if the synergistic effect of and is to be found in the formation of , catalase might have a double function (i.e., as a protective factor toward exogenous exposure of both and ).

2.2. Approach and Research Question

In this study, we develop a mathematical model of the kinetics of the key species of PTL (i.e., and ) as well as of the processes governing the interaction with a cell system, which are given in terms of the membrane diffusion rate constant and the intracellular catalase concentration. The system modeled is illustrated in Figure 1.

Figure 1

Illustration of the system representing a cell exposed to PTL.

As can be seen, the system consists of two compartments: the extracellular compartment (EC) and the intracellular compartment (IC). The two compartments are separated by the cell membrane, which some species in the system can diffuse through. Our mathematical model is explained in detail in Section 5, with references to all input data and assumptions made. Briefly, it takes into account (i) the diffusion of and from the EC (where these species are supplied by PTL) to the IC, (ii) the formation of from and (in both the IC and EC), (iii) the mitochondrial production of (in the IC), and (iv) the decomposition of and (in the IC). Furthermore, (v) since the -catalyzed consumption is considered to be the main route for -decay in biological systems (due to a high -concentration) [65,66,67,68,69], this reaction is also included.

There have already been attempts to capture the susceptibility toward exogenous of different cell lines in terms of their membrane diffusion rate constant and intracellular catalase concentration [13,14]. Two dependent variables that have been investigated recently are the intracellular steady-state concentration of and the so-called latency (which describes the reduced average reaction rates for the observed decomposition of due to the localization of encapsulated catalase in the peroxisomes). In [13], a lumped-parameter mathematical model, assuming that catalase is the major -removal enzyme, was developed and used to calculate the intracellular steady-state concentration for several cell lines. The model was calibrated to the experimental values of the measured critical parameters, and the resulting intracellular steady-state concentration was related to observed cell specific susceptibility to extracellular exposure of . The results showed that despite the fact that the experimental parameters including catalase concentration and membrane diffusion rate constant, in particular, varied significantly across cell lines, the calculated steady-state intracellular-to-extracellular ratio did not vary significantly across cell lines. In [14], it was investigated whether variations in the latency of peroxisomal catalase across cancer cell lines correlates with observed in vitro susceptibility to ascorbate at equivalent dosing of extracellular . The so-called effectiveness factor, which takes both the membrane diffusion rate and the overall reduced activity for encapsulated catalase into account, was used to quantify the effect of latency. The results suggest that latency alone is not a reliable parameter for predicting cell susceptibility to ascorbate (and hence, ).

In this study, we explore new dependent variables that could possibly explain the difference in cell susceptibility to an external addition of , with and without a simultaneous addition of , and ultimately, quantify the effect in terms of the membrane diffusion rate constant () and the intracellular catalase concentration (). Since we cannot distinguish a cancer cell from a normal cell solely by their membrane diffusion rate constant and intracellular catalase concentration, we will have to work under the notations “cancer-like cells” (i.e., systems in the higher range of membrane diffusion rate constant and the lower range of catalase concentration) and “normal-like cells” (i.e., systems in the lower range of membrane diffusion rate constant and the higher range of catalase concentration). We investigate different regimes of the supplied extracellular - and concentrations according to experimental observations of the regimes of selective/non-selective and synergistic/non-synergistic anti-cancer effect of PTL [25,26]. The dependent variables that we investigate are:

The temporal maximum ofandin the IC. As opposed to the steady-state value of the intracellular concentration, the temporal maximum can be expected to be dependent on both and . These dependent variables may be related to the maximal intracellular oxidative power of the extracellularly added (and ), and thus it would be of interest to study whether a certain extracellularly added concentration of (and ) would result in a higher oxidative power in a more cancer-like cell than in a more normal-like cell.

The system response time (i.e., the time out of equilibrium) with respect toin the IC. In order to achieve tumor progression, it is essential for cancer cells to optimize their RONS concentration and maintain the RONS equilibrium. For our mathematical model, this is translated into the question: does a more cancer-like cell have a longer response time compared to a more normal-like cell?

The “load” of intracellularand(i.e., the time integral ofandin the IC). As the temporal maximum of and cannot capture any information about the total “load” of and (i.e., how much the intracellular and is increased over a period of time), it could be of interest to study such a dependent variable as a complement. The load can be seen as a measure that combines the temporal maximum concentration and the system response time. Another possible way to define the load of intracellular would be to only consider the concentration of over a “baseline”. Here, the steady-state intracellular , before the perturbation of an addition of extracellular and at the upper limit of , is used as the baseline.

The inverse of the average and maximal rate of extracellularconsumption. Since the cell susceptibility of CAP and PTL has been found to be inversely proportional to the (extracellular) consumption rate of it is of interest to explore a dependent variable quantifying the system susceptibility in terms of the consumption. We investigate two such candidates where one was defined in terms of the inverse of the average consumption rate, and the other one in terms of the inverse of the maximal consumption rate.

For all proposed dependent variables, we will analyze the dependence on and and whether a more cancer-like cell is associated with a higher “response” than a more normal-like cell. Our main research question is thus: Can the difference in cell susceptibility toward PTL be understood, and even quantified, by one of these dependent variables?

To the best of our knowledge, this is the first study of its kind, and our aim is to take some initial steps in the direction of an increased understanding of the mechanisms underlying the selective and synergistic anti-cancer effect of PTL, and ultimately, be able to predict the response of different cells.

3. Results

As introduced in Section 2.2, in order to try to understand the combined role of the membrane diffusion rate constant and the intracellular catalase concentration in determining the susceptibility of cells toward exogenous , we have to go beyond the steady-state value of the intracellular concentration [13] (as well as latency [14]) and examine dependent variables that take the system’s temporal response of a perturbation in the EC into account. To be able to present the results in a more compact manner, the variables not yet introduced but of importance, and their denotations, are presented in Table 1. Details about the independent and dependent variables can be found in Section 5, where the mathematical model is presented. Likewise, details about the numerical calculations such as the values of the independent variables and parameters used in the model, can be found in Section 6.

Table 1

Denotations of variables used in the results analysis.

Variable	Meaning
[H_2O_2]0EC	Initial [H_2O_2] in the EC
[NO2−]0EC	Initial [NO2−] in the EC
[H_2O_2]^IC	[H_2O_2] in the IC
c_1,max	Temporal maximum of [H_2O_2] in the IC

For the analysis and interpretation of the results, we mainly consider three important features of the dependent variable of interest:

Does it account for selectivity with respect to different regimes of ?

Does it account for a synergistic effect when is added to the system?

Does it represent a feasible measure to quantify the susceptibility to exogenous of a cell system in terms of and ?

To qualify as a “measure” (i.e., as a quantification of the susceptibility in terms of and ), the dependent variable should be associated with a higher value for cells with a higher susceptibility and a lower value for cells with a lower susceptibility. Thus, in accordance with experimental observations, a feasible measure should result in a higher value for more cancer-like cells than for more normal-like cells, at least in the expected regime of selectivity (that is, for in the μM-range [26]). However, it should be noted that in this study, we do not follow strict mathematical criteria for a function to be categorized as a measure.

3.1. The Temporal Maximum of the Intracellular Hydrogen Peroxide Concentration: A Possible Measure of the Cell Susceptibility to Exogenous Hydrogen Peroxide

Our calculation results suggest that the temporal maximum of (i.e., ) is the dependent variable of major interest in terms of our requirements. Therefore, we focus our analysis on this variable. The results of the other dependent variables are presented in Appendix A.

Figure 2 shows as a function of and for with and without . The same results, but for , are shown in Figure 3.

Figure 2

The dependent variable (i.e., the temporal maximum of in the IC) as a function of and when . (a) and (b).

Figure 3

The dependent variable (i.e., the temporal maximum of in the IC) as a function of and when . (a) and (b).

When comparing the result for the different -regimes for (see Figure 2a and Figure 3a), we see that is also in different concentration regimes, which is logical. Indeed, for , whereas for , . Thus, by assuming that there exists a threshold value for which all types of cells undergo cell death, selectivity could be accounted for. However, there is no obvious synergetic effect; when comparing Figure 2a,b, is almost identical. Thus, the addition of does not change significantly. The same is true for Figure 3a,b.

For shows an increased -dependence with increasing . The lowest value of is for the lowest values of and highest values of , as would be expected for a dependent variable that would qualify as a measure of the cell susceptibility in terms of and . In addition, the highest value of is associated with the lowest value of . However, in this regime, the dependence on is insignificant. Here, in contrast, there exists a significant -dependence and by changing the scale on the -axis to a log-scale (see Figure 4), we see that there are two distinct regimes with a clear shift from one regime to another at about . The regimes of and with the most profound difference between the value of is between cells with (for all ) and cells with the lowest possible and highest possible . Thus, (i.e., the temporal maximum of ) is associated with a higher value for cancer-like cells than for normal-like cells. Indeed, is about four times greater for the most susceptible cells compared to the most resistant cells.

Figure 4

The dependent variable (i.e., the temporal maximum of in the IC) as a function of and when and .

In summary, does capture the dependence of and in a manner that is consistent with experimental observations and could thus represent a feasible measure to quantify the susceptibility of different cells in terms of their membrane diffusion rate constant and intracellular catalase concentration. However, in our model, it cannot yet account for the synergistic effect when is added.

To the best of our knowledge, there are not yet any experimental results to support our findings. We hope that our theoretical work will inspire future experimental studies. In the next section (Section 3.2), we discuss possible opportunities to experimentally quantify .

3.2. Physical Interpretation and the Use of the Temporal Maximum of the Intracellular Hydrogen Peroxide Concentration as an Experimental Probe

Our model, with all the equations, is explained in detail in Section 5. Here we use the equations to better understand how we can use as a measure to quantify the response of different cells. In order to analyze and write the equations in a more compact manner, we first introduce some short notations as well as some new notations (see Table 2).

Table 2

Denotations of variables used in the results analysis.

Variable	Meaning
c1IC	[H_2O_2]^IC
c1EC	[H_2O_2]^EC
c_2	[CATFe^III]
c_3	[CATFe^IVO^•+]
k_P	Rate of mitochondrial H_2O_2 production
k_1	Rate of H_2O_2 consumption by [CATFe^III]
k_2	Rate of H_2O_2 consumption by [CATFe^IVO^•+]

The temporal maximum of (i.e., ) occurs when the production and consumption of intracellular are equal and (as opposed to the steady-state value of , which is governed by the same rate equation but for which ). If we exclude in Equation (10) (see Section 5.2.2) the term representing the formation of from and (since it is much smaller than the other terms), we have Thus, Here, we can furthermore use the constraint since the total catalase concentration will be constant. By noting that (see Section 6.2.1), we can use the approximate expression

From our numerical calculations, we know that for low , is independent on whereas for high , is highly dependent on Furthermore (for ), for low , is highly dependent on , whereas for high , is independent on . The question is whether this behavior can be understood.

In the analysis, we first note that the implicit importance of and in determining the value of at the time of , and thus , is hidden. The dependence on originates from the fact that in our model, is determined by . Equation (3) (in Section 5.2.1) can be approximated as: Hence, the initial rate, or driving force, of -consumption in the EC will crucially depend on , and thus, . In fact, for , whereas for (see Equation (18), Section 6). It means that the initial driving force is about ten times higher in the latter case compared to the former. This could explain why is seemingly independent on at low values of ; if for all values of will be the dominant factor of the build up of in the IC. For higher values of , it seems like somewhen at , becomes increasingly dependent on . It is thus reasonable to believe that the term is becoming increasingly dominant and that the larger the value of , the larger the value of . For a fixed value of , will thus increase with increasing values of .

In summary, this means that the lower the value of , the less important is the value of , and the other way around. Thus, the susceptibility (toward exogenous ) of cancer-like cells is not much influenced by the membrane diffusion rate constant and this is due to their much higher level of intracellular prior to the perturbation by the addition of exogenous . Normal-like cells, on the other hand, are more sensitive to the value of the membrane diffusion rate constant, since the difference in concentration between the intracellular and extracellular will be much larger.

Another aspect of is whether it could provide an opportunity to extract information about different cell lines in terms of their membrane diffusion rate constant and intracellular catalase concentration. By measuring and the corresponding for different it could be possible to roughly quantify and . There are many experimental techniques for the detection and quantification of the concentration in vitro and in vivo. The intracellular concentration has been detected and measured by a chemoselective fluorescent naphthylimide peroxide probe [70], by a genetically encoded red fluorescent sensor [71], and by fluorescent reporter proteins [72]. Thus, even if does not represent a feasible measure of the cell susceptibility in terms of and , it could still possibly be used to gain more knowledge about the correlation between and and cell susceptibility toward exogenous and PTLs.

4. Discussion

In this study, we use a theoretical approach to increase the knowledge about possible underlying causes of the anti-cancer effect of PTL. Although the model is fairly simple, it does include the major pathways for species production and consumption relevant for such a cell system. It also puts emphasis on two important features (i.e., the membrane diffusion rate constant and the intracellular catalase concentration), possibly explaining the different cell responses and cell susceptibility toward PTL when comparing normal cells to cancer cells, but also when comparing resistant vs. sensitive cancer cells. Nevertheless, it is important to keep in mind that in our model, different cells are only defined in terms of these two features, which are independent variables in our analysis, whereas in reality, there are countless of other features characteristic for different types of cells that could play an important role in the context of the anti-cancer effect of PTL. Here, we merely analyze the immediate cell response determined by the scavenging system active at high concentrations of . However, we do believe that our results contribute to a better understanding of some mechanisms probably underlying the anti-cancer effect of PTL. It brings novelty to the field of plasma oncology, and more broadly to the field of redox biology by using a theoretical approach and by proposing new ways to quantify the selective and synergistic anti-cancer effect of PTL in terms of inherent factors of cells. Here, we discuss each of our main findings and their potential implications. We also highlight what we believe are the most important limitations of the model.

As opposed to the steady-state intracellular concentration of , which has been evaluated in previous studies [13], our results suggest that the temporal maximal concentration of intracellular could be a measure feasible to quantify the cell susceptibility toward exogenous in terms of the membrane diffusion rate constant and the intracellular catalase concentration. This result furthermore enables us to speculate whether the mode of action of is as a signaling molecule rather than as a toxic substance causing necrosis. It is known that the intracellular concentration of a signaling molecule rises and falls within a short period. Indeed, whether a signaling molecule is effective or not is determined by how rapidly it is produced, how rapidly it is removed, and the concentration it must reach to alter the activity of its target effector. Of particular relevance in our context is that several reports have demonstrated that the rate of generation and its concentration as a function of time play a key role in determining target cell damage or destruction [73,74,75]. RONS are regulators of signaling pathways such as the extracellular signal-regulated kinase (ERK) mitogen activated protein kinase (MAPK) pathway, which is important for cell proliferation, and a number of studies have demonstrated the ability of exogenous oxidants to activate the ERK MAPK pathway [76,77,78,79,80]. As in the general case, the duration and intensity of the ERK MAPK signal determine the outcome of the cellular response; there is a connection between the levels of ROS in a cell and the levels of MAPK signaling. In particular, MAPKs are activated in response to [81,82,83].

Based on our modeling results (presented in the Appendix A; i.e., Figure A1 and Figure A9), we do not think that the formation of itself plays a major role in the explanation of the synergistic effect of and . This is because although the overall intracellular concentration of is increased with about one order of magnitude when is added to the system, the dependence on the membrane diffusion rate constant is such that cells with a higher value of the membrane diffusion rate constant (i.e., cancer-like cells) are associated with a lower maximal intracellular concentration than more normal-like cells (i.e., cells with a lower value of the membrane diffusion rate constant). In addition, the load of is independent of the membrane diffusion rate constant. However, an important aspect to keep in mind regarding our results for and the choice to include -catalyzed consumption of in our model, is that redirects much of the produced in vivo toward radical mechanisms [65]. Indeed, many of the reactions of in vivo are more likely to be mediated by reactive intermediates derived from the reaction of with than by itself [84,85]. Thus, if the production of such reactive intermediates were to be monitored instead of , our results might be different. In this context, especially the formation of should be considered; a fraction (about 30%) of the formed will produce cage-escaped and radicals according to [86,87,88] where and [89]. A possibly important target in the context of our study is catalase; catalase is so far the best known protein target for and the rate constant of the reaction of bovine liver catalase with is at [90]. Since the temporal maximum of intracellular (i.e., ) is inversely dependent on the catalase concentration (i.e., ) with an increasingly steeper incline for lower catalase concentrations in the regime (see Figure 2), cancer-like cells would be more vulnerable to a decrease in the catalase concentration than normal-like cells, which are associated with higher values of . Thus, including these reaction pathways may possibly also make the dependent variable able to account for the synergetic effect of . Such an extension of our model was out of the scope for this study, but would be highly interesting in a future model development.

Figure A1

The dependent variable (i.e., the temporal maximum of in the IC) as a function of and when . (a) and (b).

Figure A9

The dependent variable (i.e., the load of in the IC) as a function of and when . (a) and (b).

In experiments, the consumption rate of extracellular has been found to inversely correlate with the susceptibility of cancer cells toward exogenous [32]. Thus, cancer cell lines with a high consumption rate were less susceptible. Our results cannot yet account for this correlation; when cells are defined in terms of their membrane diffusion rate constant and their intracellular catalase concentration, susceptibility in terms of the inverse of the extracellular consumption rate is not consistent with the experimental observations of cancer cells having a higher membrane diffusion rate constant and a lower catalase concentration (see Section 2.1). The fact that our model does not reproduce these patterns leaves an open question of how to construct a dependent variable in terms of the inverse of the extracellular consumption rate so that it corresponds to the experimental correlation.

The fact that our mathematical model, as well as our criteria for a dependent variable to represent a feasible measure of the cell susceptibility, does not select the system response time as a good candidate does not necessarily indicate that this variable, in general, cannot capture cell susceptibility toward exogenous . Indeed, in our definition of this dependent variable, we assume a tolerance of a 10% increase of the intracellular steady-state concentration, and a different assumption of the tolerance might give a different result.

Except for the limitations of the model already mentioned in this discussion, some other model assumptions could hamper a realistic representation of a cell system in interaction with PTL. One such limitation is that in our model, the rate of mitochondrial production is constant. Although it can be argued that this assumption is a valid starting point, in a model development, it could be important to modify this aspect to represent a cancer cell in a more realistic manner. Indeed, it has been shown that in some cancer cells, the mitochondrial respiration is decreased (in favor of aerobic glycolysis), and moreover, this shift seems to be a dynamic process (see [91] and references therein). We believe that future models could benefit from trying to take such variation of the rate of mitochondrial production into account, but this was out of the scope for this study.

Another aspect to take into account in a more realistic model is the fact that the membrane diffusion rate constant is not a static but dynamic property. In [92], it was shown that cellular stress conditions reversibly inhibit the diffusion of (and ) of aquaporin 8. Thus, a more complex model taking the implicit time-dependence of the membrane diffusion rate (caused by the increased intracellular concentration after the addition of exogenous ) could potentially produce results different from our model.

A third aspect to be aware of is that in our model, we assume that the addition of PTL does not affect the membrane diffusion rate constants. However, a number of studies have reported an enhanced cell membrane permeability (and thus, increased membrane diffusion rate constants) after CAP/PTL treatment [93,94,95]. For the aim and approach of our study, where the membrane diffusion rate constant of the key species is varied within a range of possible values, we believe that our assumption is a valid starting point. Nevertheless, for future model extensions and developments, this aspect might be important to take into account.

Finally, it should be mentioned that the rate equations used to model the system are derived from information (collected from the literature) about rate constants and reaction orders for each reaction as they appear in the experiments. Most likely, the experimental conditions will deviate from the conditions of cells treated with PTL, which will affect the accuracy of the results produced by the model. However, for the purpose of our study, we believe that parameter values of the correct order of magnitude are sufficient at this stage.

5. Mathematical Model

Mathematical models of biological reaction networks such as the system considered in this study, can generally be divided into two categories: predictive and descriptive models. Since the experimental studies on which we build our model on are primarily in vitro studies, we construct a predictive model in this work. This means that we put together the information about each of the involved reactions (reaction orders, rate constants, etc.) as they appear in the experiments. From there, the result for a certain set of initial conditions is generated by solving the time-dependent equations of motion, representing the time evolution of the system.

In this section, we systematically present the species and reactions in the system considered (Section 5.1), how the system time evolution is modeled (Section 5.2), and we explicitly define the dependent variables that are analyzed (Section 5.3).

5.1. Species and Reactions in the System

The involved species of interest are , , , , and (see Figure 1 in Section 2). The following reactions and interactions of the species , and (native) catalase in the system are taken into account.

Decomposition of by catalase:

Generation of through reaction between and : where the equilibrium between and is described by:

Decomposition of by catalase:

-catalyzed consumption of :

The denotations of the time-dependent concentrations of the different species are shown in Table 3.

Table 3

Denotations of the time-dependent concentrations in the system.

Species	Denotation
[H_2O_2]	c_1
[CATFe^III]	c_2
[CATFe^IVO^•+]	c_3
[ONOO^−]	c_4
[NO2−]	c_5
[H^+]	c_6
[ONOOH]	c_7
[CO_2]	c_8

5.2. Modeling the System

The mathematical model considers the kinetics of the reactions in the system composed of two subsystems (EC and IC), see Figure 1 in Section 2, as well as diffusion of certain species between the two subsystems. The equation governing the kinetics of each species is given by the sum of the reaction rates (describing the rate of production and consumption of species ), and (in the case of species 1, 4, 5, and 7), the diffusion rate through the cell membrane, from the EC to the IC,

Equation (1) represents the resulting rate equation, derived from the rate constants and reaction orders for each reaction as they appear in the experiments. Equation (2) describes the rate of membrane diffusion of species according to Fick’s law of diffusion with a linear concentration gradient over the cell membrane. Here, is the rate of species exchange through the membrane. We denote this as “membrane diffusion rate constant”. More information about the derivation of Equation (2) can be found in Appendix B.

Explicitly, our mathematical model is used to analyze the behavior of a dependent variable , where denotes the set of independent variables that are varied in the system. The independent variables in our model are:

The membrane diffusion rate constant through the cell membrane () and

The initial intracellular catalase concentration ().

The species and are assumed to be present in equal initial concentrations in both the EC and the IC (thus, and ). Since we do not explicitly study the kinetics of these species, we make such an assumption to reduce the complexity of the model.

Detailed information about the mathematical model is presented in the following sections.

5.2.1. Mathematical Model of the Reaction Kinetics in the Extracellular Compartment

At and in certain initial concentrations ( and ) are inserted into the EC, representing treatment of the cell by PTL (as these species are the dominant RONS in PTLs), and their reactions as well as diffusion through the membrane into the IC is monitored. The reaction network and resulting set of differential equations are given below.

Reaction network

Differential equations

5.2.2. Mathematical Model of the Reaction Kinetics in the Intracellular Compartment

At the concentration of is at a certain steady-state value () because it is continuously produced by the mitochondria at the rate and decomposed by catalase, which exists in the IC, and is modeled as free in the solution. (More information can be found in Section 6.2.3). It is assumed that at the total amount of catalase exists as (i.e., . The reaction network and resulting set of differential equations are given below.

Reaction network

Differential equations

Equations (3)–(8) and (10)–(17) are solved numerically. The details about the numerical calculations can be found in Section 6.3.

5.3. Dependent Variables

In the following sections, we explicitly define the dependent variables analyzed in this study.

5.3.1. Temporal Maximum of the Intracellular Hydrogen Peroxide and Peroxynitrite Concentration

The dependent variables and are defined as and

5.3.2. System Response Time of Intracellular Hydrogen Peroxide

Assuming that the system has a tolerance of an increase of 10% of the baseline concentration (see Section 3.2), the dependent variable can be formulated

5.3.3. Load of Intracellular Hydrogen Peroxide and Peroxynitrite

The simplest way of creating a quantitative measure of the “load” of intracellular and is to use the time-integral over the whole time regime () as the dependent variable, in other words, and

For the “load” over the baseline concentration of intracellular , if we denote this baseline constant the dependent variable is defined as

5.3.4. Rate of Extracellular Hydrogen Peroxide Consumption

Here, we first define the dependent variable as The average extracellular consumption rate of is then defined as and the maximal extracellular consumption rate of as

In order to create a potential measure (i.e., a dependent variable where a more cancer-like cell is associated with a higher susceptibility), we use the variables and in our calculations.

6. Numerical Calculations

6.1. Independent Variables

The two independent variables in the system are and (i.e., the diffusion rate constant of through the cell membrane from the EC to the IC, and the initial catalase concentration in the IC). Furthermore, we use four different combinations of and in our calculations. The motivation and details of these variables are given in the following sections, and have also been introduced in Section 2.

6.1.1. Membrane Diffusion Rate Constant of Hydrogen Peroxide

In [96], the diffusion rate constant for crossing lipid bilayers was found to be . Due to the chemical similarities between and , we use this value as a reference value for and we vary within the range .

6.1.2. Initial Catalase Concentration in the Intracellular Compartment

The intracellular concentration of catalase is calculated from two different premises (see Appendix C). Considering the rough estimates in both approaches, it seems reasonable to use an effective catalase concentration in the range of in our calculations. As a reference value, catalase concentration in human blood cells is about [97,98].

6.1.3. Initial Hydrogen Peroxide and Peroxynitrite Concentration in the Extracellular Compartment

Several publications have shown that and are formed at concentrations ranging from μM to mM in plasma-treated liquids (PTLs) [99,100,101,102]. In this study, we use the initial conditions for and shown in Table 4. The different regimes of these four combinations are specified in the last column. We assume that the selectivity is related to the concentration of extracellular (i.e., selective cancer killing only occurs at low concentrations (order of [26]), while higher concentrations (e.g., order of ) kill both cancer and normal cells) [25]. Based on this assumption, we want to compare the dependent variables for the selective versus non-selective regime. For both regimes (selective versus non-selective), we furthermore want to investigate whether a synergistic effect can be found (i.e., if the values of the dependent variables are enhanced when and are added together) [25,26].

Table 4

Initial concentrations of and in the extracellular compartment.

[H_2O_2]^EC (M)	[NO2−]^EC (M)	Regime
10^−3	10^−3	Non-selective, synergistic
10^−3	0	Non-selective, non-synergistic
10^−6	10^−3	Selective, synergistic
10^−6	0	Selective, non-synergistic

6.2. Parameter Values

6.2.1. Reaction Rate Constants

The used rate constants are summarized in Table 5, along with the references where the data is adopted from, and some remarks about the conditions for which these values were reported.

Table 5

Reaction rate constants.

Rate Constant	Parameter Value	Reference	Remark
k_1	1.7×10^7 M^−1s^−1	[103]	Mammalian catalases
k_2	2.6×10^7 M^−1s^−1	[103]	Mammalian catalases
k_3	1.1×10^3 M^−2s^−1	[37]	At pH=3.3 and T=25 °C
k_4	K_ak_−4=10^−pK_ak_−4	[104,105]	The pK_a-value at T=25 °C is 6.5–6.8
k_−4	∼10^10 M^−1s^−1	[106]
k_5	1.7×10^6 M^−1s^−1	[64]	At pH=7.1 and T=25 °C
k_6	5.8×10^4 M^−1s^−1	[84]	At T=37 °C, pH-independent

6.2.2. Membrane Diffusion Rate Constants

, when protonated (i.e., as ), is reported to diffuse easily across biological membranes [107]. When not protonated, anionic channels have been shown to be permeable to [108]. It has furthermore been established that is able to penetrate cell membranes [96,109]. In [96], using model phospholipid vesicular systems, it was demonstrated that freely crosses phospholipid membranes. The diffusion rate constant for crossing lipid bilayers was found to be . Due to the acid–base equilibrium between and its conjugated acid , this is likely an average value for and . Thus, . Since is an anion as well as similar in size, we assume the same value (i.e., ). The used diffusion rate constants are summarized in Table 6. Note that we do not consider the potential effect of the PTL on the membrane diffusion rate constants in our model.

Table 6

Membrane diffusion rate constants.

Rate Constant	Parameter Value	Reference	Remark
k_D,4	320 s^−1	[96]
k_D,5	320 s^−1		Assigned
k_D,7	320 s^−1	[96]

6.2.3. Initial Concentrations

The used initial concentrations in the EC and IC are summarized in Table 7 and Table 8, respectively.

Table 7

Initial concentrations of the species in the extracellular compartment.

Species	Initial Concentration (M)	Reference	Remark
H_2O_2	Varied
ONOO^−	0		Assigned
NO2−	Varied
H^+	10^−7		Assigned
ONOOH	0		Assigned
CO_2	10^−3	[110]

Table 8

Initial concentrations of the species in the intracellular compartment.

Species	Initial Concentration (M)	Reference	Remark
H_2O_2	Varied
CATFe^III	Varied
CATFe^IVO^•+	0		Assigned
ONOO^−	0		Assigned
NO2−	10^−4	[111,112,113]	See Appendix D
H^+	10^−7		Assigned
ONOOH	0		Assigned
CO_2	10^−3	[110]	See Appendix D

The initial concentration of intracellular (i.e., ) is varied with the initial concentration of catalase in order to achieve the correct steady-state for each The generation from mitochondria is in the range of [114], which corresponded to [115].

Thus, from Equation (10), at (and thus, the term vanishes), Assuming that at steady-state, Hence, in our model.

6.3. Software and Details about the Calculations

The numerical calculations are performed in MATLAB. Due to significant differences in time scales, we use the solver ode23s to solve the set of rate equations.

The simulations are performed at time-scales covering the transient of the system’s response. For the calculations, we use the time intervals and time steps shown in Table 9. We start with a very short time-step in the first 10 ms, which is then enlarged by a factor of 100 until 1 s, and again by a factor of 100 until the final time of 100 s.

Table 9

Time intervals and time steps.

Time	Value (s)	Time Step	Value (s)
t_1	10^−2	dt_1	10^−7
t_2	1	dt_2	10^−5
t_f	10^2	dt_3	10^−3

We furthermore vary the independent variables according to Table 10.

Table 10

Minimal and maximal values as well as number of steps of independent variables.

Independent Variable	Minimal Value	Maximal Value	Number of Steps
[CATFe^III]	10^−8 M	10^−5 M	100
k_D,1	100 s^−1	2000 s^−1	100

7. Conclusions

With this study, we aim to gain insights about the possible mechanisms underlying the anti-cancer effect of plasma-treated liquids (PTLs). In particular, we are interested in whether cell susceptibility toward PTL can be quantified in terms of cell-specific features, how selectivity arises, and why combined with (as in PTL) offers a synergistic and thus enhanced anti-cancer effect compared with only. By developing a mathematical model describing the kinetics of the species in PTL-treated cells, we analyze four different dependent variables as a function of the membrane diffusion rate constant and the intracellular catalase concentration. Ultimately, one or more of these dependent variables could be used to quantify selective and synergistic effects of PTLs for different types of cells. In accordance with experimental observations, cancer cells are supposed to be associated with a higher membrane diffusion rate constant and a lower intracellular catalase concentration compared to normal cells, and we use this knowledge in the evaluation of our proposed dependent variables.

The model is built up ab initio based on the species, reactions, and processes of major importance in the context of cell susceptibility toward PTL, and parameter values such as rate constants are extracted from the literature. Thus, the model itself summarizes the current state of knowledge on the matter in a compact and descriptive manner. This type of mathematical modeling to gain insight into the underlying mechanisms of the anti-cancer effect of PTL is novel and this study is the first of its kind in the field of plasma oncology. Furthermore, we propose new ways to quantify the selective and synergistic anti-cancer effect of PTL in terms of inherent cell features, which is also an innovative approach in the ongoing research on the mode of action of PTL.

We find that the temporal maximal intracellular concentration shows a dependency of the membrane diffusion rate constant and the intracellular catalase concentration, so that it could possibly be used to quantify the anti-cancer effect of exogenous , but it does not account for the synergistic effect of and in PTL. However, by including the reactions where is produced in the catalyzed consumption of , and the interaction between and catalase, the dependent variable could possibly also be able to account for the synergetic effect of .

We believe that our model is an important step to unveil the underlying mechanisms of the anti-cancer effect of CAP and PTLs, but more efforts are needed in order to understand the full picture of the causes and action. Here, both positive and negative results are important to share, in order to increase our collective knowledge of which clues may lead us forward in our search, and which clues we can leave behind, at least for now. Theoretical and experimental approaches to investigate possible key features of cells and their interaction with CAP and PTLs play complementary roles in our aim to push the limit of knowledge further. We hope, and believe, that our study contributes to the quest to quantify selective and synergistic effects of plasma for cancer treatment.