Enhanced release of primary signals may render intercellular signalling ineffective due to spatial aspects

Detailed mechanistic modelling has been performed of the intercellular signalling cascade between precancerous cells and their normal neighbours that leads to a selective removal of the precancerous cells by apoptosis. Two interconnected signalling pathways that were identified experimentally have been modelled, explicitly accounting for temporal and spatial effects. The model predicts highly non-linear behaviour of the signalling. Importantly, under certain conditions, enhanced release of primary signals by precancerous cells renders the signalling ineffective. This counter-intuitive behaviour arises due to spatial aspects of the underlying signalling scheme: Increased primary signalling by precancerous cells does, upon reaction with factors derived from normal cells, produce higher yields of apoptosis-triggering molecules. However, the apoptosis-triggering signals are formed farther from the precancerous cells, so that these are attacked less efficiently. Spatial effects thus may represent a novel analogue of negative feedback mechanisms.


Hydroxyl radicals
Hydroxyl radicals ( • OH) are extremely short-lived (highly reactive) species; therefore, their diffusion has been neglected here. From the reaction-diffusion equations (Eq.5 in the main text), their quasi-steady-state concentration profile can be derived, for which [ OH • ]( , ) ≈ 0 over timescales T much longer than OH but short compared to the lifetimes of longer-lived species and to characteristic times of cellular processes such as proliferation or apoptosis: Here the first term describes the production of hydroxyl radicals upon the reaction of hypochlorous acid with superoxide (peroxidase pathway). The second term captures • OH from peroxynitrite pathway, i.e. from the decay of ONOOH, the conjugate acid to peroxynitrite anion ONOO -, with lifetime ONOO and efficiency OH~3 3 % [Lobachev and Rudakov 2006] where k LPO = 10 9 M -1 s -1 [Radi et al 1991] stands for the reaction rate constant for lipid peroxidation by • OH, n lipid ≈ 10 -15 mol/cell [Alberts et al 2002) for the amount of lipid molecules in cell membrane that could be attacked, and ( , ) for cell density at time T, at x = 0 for transformed cells and at x = L 1 ~ 1 mm for effector cells in co-culture experiments.

Hypochlorous acid
Hypochlorous acid is also relatively short-lived, quickly reacting with constituents of intercellular medium [Kundrát et al 2012, Deborde andvon Gunten 2008]. Hence its diffusion can be neglected with respect to its lifetime, and for time scales T long compared to HOCl lifetime it holds where superoxide concentration at x=0 has been taken as an upper limit on superoxide levels (see next subsection), somewhat overestimating the effect of the reaction on HOCl lifetime. The absorption of HOCl in reaction with H 2 O 2 and POD has been neglected. Assuming peroxidase be a relatively stable enzyme [Kundrát et al 2012], it is distributed almost homogeneously, [POD]( , ) ≈ [POD]( ).

Superoxide
For steady-state concentration profile of superoxide, the reaction-diffusion equation ( The boundary conditions are where O 2 −• stands for the release rate of superoxide per transformed cell (assumed to be constant) and ( ) for the density of transformed cells, which is time-dependent due to cell proliferation and apoptosis.
Neglecting the effect of reactions and assuming that the diffusion length of superoxide O 2 −• = √ O 2 −• O 2 −• is much smaller than the height of cell culture medium, O 2 −• ≪ ~ 3 mm, the steady-state concentration profile of superoxide is In particular, at x=0 where transformed cells are,

Nitric oxide
Nitric oxide is produced by transformed cells (at x=0, cell density , per-cell release rate NO • ) as well as by effector cells (at x=L 1 , cell density , per-cell release rate NO • ), so that the boundary conditions read where + and − denote right and left derivatives.
If the effect of reactions is neglected and the diffusion length NO • = √ NO • NO • is much smaller than the height of cell culture medium, NO • ≪ , we obtain Note the factor of 1/2 in the second term, which accounts for the diffusion from the source at x=L 1 into two directions.

Hydrogen peroxide
For hydrogen peroxide the quasi-steady-state approximation to the reaction-diffusion equation reads This result could be interpreted as if hydrogen peroxide were produced by dismutation of superoxide locally at x=0 only and accumulated over the period of * H 2 O 2 . Note however the non-trivial factor 1 + 2 * H 2 O 2 O 2 −• that accounts for the diffusion properties of superoxide and its dismutation product, hydrogen peroxide.

Antioxidants converting superoxide to hydrogen peroxide
The concept of superoxide lifetime used above accounts for reactions not explicitly included in the given reaction scheme that remove superoxide from the signalling. Numerous antioxidant systems, however, act via converting superoxide to hydrogen peroxide. This is the case e.g. for superoxide dismutase (SOD) or vitamin C [Gray andCarmichael 1992, Saran andBors 1994].
With such antioxidants, e.g. SOD, the time-and space-dependent concetration of superoxide is given by where O 2 −• +SOD denotes the corresponding reaction rate constant, e.g. 6.4×10 9 M -1 s -1 for Cu,Zn-SOD [Gray and Carmichael 1992]. For a spatially homogeneously distributed antioxidant, its removal of superoxide can be included in the effective lifetime of superoxide, The first term describes the formation of hydrogen peroxide from superoxide by spontaneous dismutation; this process is quadratic in superoxide concentration, and hence leads to the factor of 2 in denominator. The second term describes dismutation catalysed by SOD (or other antioxidants); this conversion is linear in superoxide concentration, and hence the factor of 2 in denominator is not present.

Mutual reactions approximated by local absorption
In addition to having assumed that the diffusion length of superoxide be much smaller than the medium height, i.e. as a product of the reaction rate constant, the concentrations of the reactants, and the effective size of the reaction region, which equals O 2 −• for superoxide dismutation.
Apart from NO • autoxidation, the equations for NO • are linear, so that NO • derived from transformed and from effector cells could be treated separately. Similarly to the above-discussed case of superoxide dismutation, for the terms describing the reaction of superoxide with nitric oxide derived from transformed and effector cells, respectively, it holds (with | . | denoting absolute value of the argument): the distinct signs in denominators follow from the exponentials in the two contributions to the concentration profile for nitric oxide, Eq.8.
Thus, instead of including these reaction terms in the reaction-diffusion equations directly, they could be approximately accounted for via boundary conditions. In this local absorption approximation, the quasi-steadystate concentration profiles of superoxide and nitric oxide are thus given as solutions to the following set of differential equations and boundary conditions: Equations 15-18 represent a set of coupled equations: Superoxide concentrations at x=0 and x=L r obtained from Eq.15 affect the absorption terms γ 0 and γ r in boundary conditions of Eqs. 16-17 for nitric oxide, and vice versa, nitric oxide levels influence superoxide concentrations through the absorption terms β 0 and β r . The two species are also linked by Eq.18 that defines the effective reaction point L r .
Two methods of solving this set of equations are described below. The first method is a relatively crude approximation, neglecting the competition for superoxide between the three signalling modes (peroxidase pathway and the peroxynitrite pathway with NO • derived from transformed and effector cells, respectively).
However, it provides analytical formulas for quasi-steady-state concentration profiles of the signalling species and for the efficiencies of the distinct signalling modes relevant for intercellular induction of apoptosis. The second method accounts for the interplay of the pathways and modes via an iterative procedure.

Analytical formulas for quasi-steady-state concentration profiles and efficiencies of signalling modes
For constant absorption terms β 0 , β r , γ 0 and γ r , analytical solutions can be found to Eqs.15a-d, 16a-c and 17a-e.
Comparing these formulas with Eqs.6 and 8 reveals that the absorption terms at x=0 and x=L r act as if they reduced the release rates of superoxide by (1+β 0 ) and (1+β r /2), respectively; similarly for nitric oxide.

Hypochlorous acid pathway
For the HOCl pathway, neglecting in Eq.15 the reactions with NO • , i.e. the terms containing [NO • ] (0, ) and [NO • ] ( , ), the absorption term β 0 contains only the dismutation term, Putting this into Eq.19a, we obtain a quadratic equation Using the approximation √1 + − 1 ≈ min (√ , 2 ), illustrated in Fig.S1 (with min( . ) standing for the minimum of the arguments), we may also write Combining Eqs.1, 2, 3, 9b, 11 and 20, we obtain the following analytical formulas for the concentration of hypochlorous acid: and for the efficiency of the HOCl pathway in terms of the levels of hydroxyl radicals attacking transformed cells:

Peroxynitrite pathway
For the ONOO pathway, the approximation of the reaction between superoxide and nitric oxide by local absorption, at x=0 or at x=L r for nitric oxide derived from transformed and effector cells, respectively, also means that peroxynitrite is modelled as if produced locally at these points. The corresponding fluxes of peroxynitrite at x=0 or x=L r read where we have assumed, similarly to the cases of superoxide and nitric oxide, that the diffusion length of peroxynitrite is negligible with respect to the medium height, ≫ ONOO .
Similarly, for the contribution to the peroxynitrite pathway from nitric oxide released by transformed cells, To be able to evaluate peroxynitrite levels in Eqs.23a-b, we thus need only the concentrations of superoxide and nitric oxide, at x=0 or at x=L r for nitric oxide derived from transformed and effector cells, respectively. These are given by Eqs. 19a-19e. Unfortunately, even if the consumption of superoxide by its dismutation were neglected, Eqs. 19a-19e could not be solved analytically in their full form. Yet, if each pathway is considered separately, i.e. if their interplay through competing for superoxide is neglected, analytical results could be derived:

Autocrine mode of peroxynitrite pathway
For the autocrine mode of the ONOO pathway, i.e. both superoxide and nitric oxide derived from transformed cells, neglecting the interplay between the pathways/modes means putting β r =γ r =0 and neglecting the dismutation term in β 0 (Eq.15e). Combining Eqs.19a, 15e and 19c, we obtain a quadratic equation for superoxide levels at the transformed cells, with a unique positive solution Similarly, from Eqs.19c, 16d and 19a, the concentration of nitric oxide reads The concentration of peroxynitrite formed is given then by combining Eqs.24a-b with Eqs.22a and 23b. Applying (24c) this approximation works in particular due to the large value of the reaction rate constant O 2 −• +NO • , as superoxide and nitric oxide react virtually whenever they encounter each other (almost diffusion-limited reaction). The amount of peroxynitrite formed is thus given directly by the production rates of superoxide or nitric oxide by transformed cells, whichever is smaller.

Inter-culture mode of peroxynitrite pathway
Neglecting β 0 in Eq.19b and requiring the amounts of superoxide and nitric oxide be equal at x=L r (which is the definition of L r , Eq.19e), Eq.19b yields a quadratic equation for [O 2 −• ]( , ), with a single positive solution, Similarly, for the concentration of nitric oxide released from effector cells that is present at x=L r one can derive where ln( . ) denotes natural logarithm.
Combining this relation with Eqs.23, 25 and 26, after simplifying we obtain This formula captures the case when the concentrations of superoxide and nitric oxide equal at some point between the cultures of transformed and effector cells, i.e. the case of 0 ≤ ≤ 1 . If less superoxide is produced than the amount of nitric oxide released by effector cells that diffuses to transformed cells (so that the 'meeting point' would be < 0) or if less nitric oxide is produced than superoxide diffusing to effector cells (so that > 1 ), Eq. 28' cannot be derived in this form as the underlying assumption that √1 + − 1 ≈ √ does not hold in Eqs.25-26. Eq.28' would provide unrealistically high results in such cases. A simple solution to this problem is considering, similarly to the case of peroxynitrite from nitric oxide released by transformed cells discussed above, that not more peroxynitrite can be formed than the levels of superoxide and nitric oxide present at transformed cells: Note that the first two terms are directly proportional to the densities of transformed and effector cells, respectively. On the other hand, the third term depends on cell densities (and per-cell release rates of superoxide and nitric oxide) in a highly non-trivial manner: Typically nitric oxide is significantly more stable than peroxynitrite, so that NO • > ONOO , and hence the first exponent is negative. This means that the yields of peroxynitrite at x=0 decrease with increasing density of transformed cells or per-cell release of superoxide. The reason for this inversed behaviour is the following: Keeping the production of nitric oxide by effector cells at L 1 constant and increasing the production of superoxide by transformed cells at x=0, the effective reaction point gets farer away from the transformed cells, i.e. L r increases. More superoxide and nitric oxide are present and hence more peroxynitrite is formed at this new reaction point. However, this effect is outweighed by the increased distance peroxynitrite has to diffuse to reach transformed cells, so that its levels there actually decrease.

Analytical approach: Summary
Taken together, the analytical approach consists in calculating peroxynitrite levels by Eqs.28 and 24c for nitric oxide derived from effector cells and transformed cells, respectively, the concentrations of hypochlorous acid by Eq.21a, and finally by converting these to signalling efficiency in terms of the yields of hydroxyl radicals by Eq.1.
The analytical formulas derived above are based on three assumptions: (1) Reactions can be approximated by local absorption; (2) Interplay between the pathways and modes can be neglected; Diffusion lengths of all species are much smaller than the height L of intercellular medium (the size of the region of interest), distance L 1 between cell cultures, and the distance L r between the transformed population and the effective source of peroxynitrite.
The assumption (3) has been introduced only to keep the analytical formulas relatively simple; it is equivalent to assuming that the species be relatively short-lived. Full formulas could be derived without this assumption but are not reported here for the sake of simplicity. On the other hand, abandoning assumption (2) of independent pathways leads to sets of coupled equations that cannot be solved analytically. The assumption (1) is crucial for the modelling presented here, although e.g. a reaction-diffusion equation for superoxide with its dismutation only but without any further reactions could be solved in terms of special functions. Note that the assumption (1) is almost perfectly fulfilled due to the high reaction rate constants in the given signalling cascade, in particular the one for the reaction of superoxide with nitric oxide. The approximation by local absorption, assumption (1), works however surprisingly well (cf. the results presented in the main text) also for superoxide dismutation, whose reaction rate constant is 4 orders of magnitude lower; the short lifetime of superoxide assures this.

Iterative method
A significant improvement to the above discussed analytical formulas can be obtained if the defining Eqs.15-18 are not solved analytically with numerous approximations as discussed above but iteratively using the perturbation theory: The In an analogous way the iterative procedure is followed to higher-order terms. In this paper the iterative procedure has been stopped when the concentrations of superoxide and nitric oxide at the effective reaction point L r differed by less than 0.01 %. If this criterion was not met, indicating a poorly convergent procedure, the procedure was limited to 500 iterations in order to restrict the needed computation time.
Contrary to the analytical formulas, in the iterative procedure the assumption (3)  Obviously, not neglecting these factors is especially important when modelling the reaction system for small distances between the populations of transformed and effector cells and/or for low amounts of medium.

Approximating mutual reactions by local absorption
Mutual reactions of signalling species can be approximated by local absorption. This is illustrated in Fig.S2 for the reaction between superoxide and nitric oxide for three values of superoxide release rate per transformed cell, namely 10 -18 mol/s, 10 -16 mol/s (the standard value from Tab.3), and 10 -14 mol/s. Quasi-steady-state concentration profiles obtained by detailed numerical simulations are depicted by histograms. Note that nitric oxide (green lines) has been generated from normal cells (at x=1 mm) as well as from transformed cells (at x=0), whereas only transformed cells produce superoxide (red lines). The iterative procedure reasonably approximates the concentration profile of normal cell-derived nitric oxide from its source at x=1 mm (and, similarly, of superoxide from its source at x=0) down to the effective reaction point L r (about 0.04, 0.35, and 0.5 mm, respectively, for the three release rates of superoxide). Outside these regions, the iteratively calculated concentration profiles underestimate the effect of reaction; they decrease with the same linear coefficient (extinction coefficient) as if no reaction took place. However, as only species concentrations at x=L r are needed for the formula that provides the reaction flux in this approximation, the predicted yields and spatially-dependent concentrations of peroxynitrite (blue lines) are reproduced satisfactorily. Even an overestimation of the effective reaction point for the highest release rate of superoxide (Panel C) does not affect peroxynitrite levels at transformed cells (x=0).
Note also that the shift of the effective reaction point L r farer away from the transformed cells with increasing superoxide release by transformed cells means that more peroxynitrite diffuses to effector cells at x=1 mm (

2.2
Roles of parameters affecting the signalling efficiency The effectiveness of IIA signalling, assessed in terms of the yields of hydroxyl radicals at transformed cells, depends sigmoidally on the release rate of NO • by effector cells (Fig.S3 Panel A). As soon as this rate exceeds about 2×10 -17 mol s -1 per cell (for other parameters taking their standard values from Tab.3), the peroxynitrite pathway with NO • derived from effector cells (green line) dominates the signalling. Analytical and iterative calculations (dashed and solid blue lines) correctly reproduce the maximal signalling efficiency effectiveness but overestimate the effectiveness of high NO • releases; this is due to having neglected the autoxidation of NO • , reactions #9-11 in Table 1, as demonstrated by the results of numerical simulations neglecting these reactions (empty squares).
Simulations with the standard parameters are insensitive to the lifetime of NO • (Fig.S3 Panel B), as this parameter does not critically influence the peroxynitrite pathway in its autocrine mode (cf. Eq.24c) and the inter-culture mode is too weak (green line). However, at higher release rates of NO • from effector cells where the peroxynitrite pathway in its inter-culture mode plays an important role, an increase in the lifetime of NO • does enhance the signalling effectiveness (simulation results shown by empty symbols); again, the iterative approach (dash-dotted blue line) overestimates the simulation results due to having neglected the autoxidation of NO • .
Similar overall sigmoid behaviour is obtained for varying the level of peroxidase present ( Fig.S3 Panel C) and the lifetime of HOCl (Fig.S3 Panel D). The outcome of the signalling is in the studied ranges only mildly affected by the lifetimes of hydrogen peroxide or peroxynitrite (Fig.S3 Panels E-F). Geometrical factors such as the distance between the cultures or the height of cell culture medium influence the signalling too ( Fig.S3 Panels G-H). The analytical formulas provide correct trends, except for the role of geometrical factors, due to having neglected in the formulas factors such as 1-exp(-L/r NO ). The iterative approach accounts for these factors, and reproduces the simulation results correctly. In fact, the results of the iterative approach are systematically slightly higher than the simulation results. This issue, however, seems to be related to a slight inaccuracy of the numerical simulations, as indicated in Fig.S3 Panel G by the results of exemplary numerical simulations using twice smaller spatial grid (red points) or time steps (yellow points). Due to their high computational expensiveness and only minor differences in the results, such simulations with finer steps have been limited to the few cases presented here.
The results of simulations and modelling on the effect of antioxidants such as superoxide dismutase (SOD) that convert superoxide into hydrogen peroxide are presented in Fig.S3 Panels I-J. In Panel I, increasing levels of SOD have been considered in addition to the standard rate of superoxide removal, described by its lifetime of 1.7 s. The simulation results (points) are nicely reproduced by the iterative approach (solid blue line), with a slight overestimation at high SOD levels; the analytical approach (dashed blue line) shows somewhat larger deviations from the numerical simulations. Note that at low and medium SOD levels the HOCl pathway (red line) dominates, while at high SOD levels the major contribution to IIA signalling comes from the peroxynitrite pathway with NO • derived from effector cells (green line). The shapes and mutual positions of these efficiency curves for individual signalling modes vary with system parameters. Enhancing the lifetimes of superoxide and nitric oxide and reducing the level of POD separates the two peaks and forms a local minimum in overall IIA signalling efficiency at SOD concentration of about 2×10 -10 M (  Table 1. In Panel B, empty symbols and dash-dotted blue line capture simulations and iterative results with NO • release rate per effector cell enhanced to 9×10 -17 mol/s. In Panel G, the cell culture distance L 1 has been varied, keeping the height of the cell culture medium at L = 3 mm. Also presented are results of simulations with twice smaller spatial grid (red symbols) or twice smaller time steps (yellow symbols). In Panel H, the height of the culture medium has been varied, keeping the inter-culture distance L 1 = L/3. In Panel I, in addition to the non-specific removal of superoxide with the standard lifetime of 1.7 s, its conversion to H 2 O 2 by SOD has been considered. In Panel J, results of exemplary calculations with the iterative approach are shown with lifetimes of superoxide and nitric oxide increased three and ten times and peroxidase level reduced ten times, other parameters unchanged.

Roles of basic cellular parameters
In Figure S4 is shown the effect of varying model parameters that describe basic cellular properties. Allowing transformed cells to reach higher densities somewhat enhances the percentage of apoptotic cells in the first maximum but leaves the long-term percentage largely unaffected ( Figure S4A). On the contrary, allowing normal cells to grow to higher densities enhances the available nitric oxide and peroxidase and makes virtually all transformed cells undergo apoptosis ( Figure S4B). Reducing the rate of removal of apoptotic cells, i.e. increasing the characteristic time t rm of this process, leaves more apoptotic cells present and, hence, increases their percentage ( Figure S4C). The rate and extent of apoptosis in the transformed population also increase if cellular repair of membrane damage gets compromised, i.e. if characteristic time for repair t rep increases ( Figure S4D). Enhancing cellular sensitivity to membrane damage, i.e. reducing the amount of damage that triggers apoptosis, shifts the onset of apoptosis towards earlier times and enlarges its overall extent ( Figure S4E). Increasing the slope of the non-linear cell response to membrane damage also affects the rate and extent of apoptosis; it may make the kinetics of apoptosis quickly oscillating ( Figure S4F). Varying the characteristic duration of apoptosis execution naturally also modulates the kinetics of apoptosis manifestation ( Figure S4G). Faster proliferation of transformed or effector cells affect apoptosis as shown in Figure S4H