### Model results for a changing environment

In this section, we take a closer look into how the mechanisms we have identified (*α*, *b*_{Apoptosis}, *c*_{Necrosis}, *m* and *n*) contribute to the final distribution of necrotic and apoptotic strategies in the system. In order to do this, we explore how nature plays the evolutionary game to arrive at optimal fractions of neutrophil population when certain parameters are varied. The motivation for this is that certain factors that affect the cost and benefit of choosing either strategy is assumed dependent on the type of environment the human body is exposed to. We have identified the following parameters *α*, *b*_{Apoptosis}, *c*_{Necrosis}, *m* and *n* as determining factors to the fitness of neutrophils. However, due to a high degree of freedom that our model exhibits, a global sensitivity analysis would not deliver meaningful results. Hence, we vary the parameters one at a time for different values of *m* and look at the effect it has on the stable distribution of strategies in the system.

We start by looking at the optimal fractions of necrotic populations (*p*) by simultaneously varying the cost of necrosis (*c*_{Necrosis}) and the amount of ITMs a single apoptotic neutrophil can neutralize (*n*) in Fig. 3 (*left column*) while fixing parameters *α*, *b*_{Apoptosis}, and *m*. Using our model, we show that the system evolves to an optimal fraction of necrotic and apoptotic neutrophils that corresponds to an increasing cost of necrosis as well as increasing capacity of ITM neutralization by apoptotic neutrophils as the insult intensifies. For the system to survive increasing threat due to ITMs, the system to increase its ITM resolution capacity by increasing *n* while also making it costly to go into necrosis. Necrosis, in the biological point of view, has detrimental effects to the system brought about by aggravating the initial level of inflammation by generating more ITMs.

Next we vary the global cost (*α*) and benefit of apoptosis (*b*_{Apoptosis}) while fixing *c*_{Necrosis}, *n*, and *m* as shown in Fig. 3 (*middle column*). Here we show that the system copes with increasing levels of insult by shifting the optimized fractions of neutrophils towards decreasing global cost and increasing the benefit of apoptosis. With all other parameters held constant, that is with a fixed ITM neutralizing capacity of both apoptotic and necrotic neutrophils coupled with also a fixed cost of necrosis, nature plays an evolutionary game where the system copes with an increasing level of insult by lowering the cost imposed by the remaining ITMs and also by increasing the benefit that can be reaped from the act of apoptosis via anti-inflammatory factors.

Finally, in Fig. 3 (*right column)*, we vary *b*_{Apoptosis} and *c*_{Necrosis} while fixing the values for *n*, *m* and *α*. Here we show that with increasing intensity of insult, the optimized fractions of neutrophils are achieved after playing the game only with increasing contributions from the cost of necrosis and benefit of apoptosis. Hence, by increasing the cost to go into necrosis coupled with increasing the anti-inflammatory effects of apoptosis, the game finds stable concentrations of neutrophil populations.

There is an apparent transition of strategy from largely necrosis to apoptosis (that is, from dominantly red to blue from left to right) as *n* (left column), *b*_{Apoptosis} (middle column) and *c*_{Necrosis} (right column) increases, thus favoring apoptosis. Hence, we interpret that increasing the cost of necrosis lowers the fitness for the necrotic strategy. Therefore, the system stabilizes into apoptosis (apparent in left and right columns). Conversely, increasing the cost of remaining ITMs (*α*) in the system shifts the equilibrium to necrosis. Note that we also see the same trends when *α* is plotted with respect to *c*_{Necrosis}, the same as what we see in the middle column. Hence, with increasing initial concentrations of ITMs, the stable evolutionary strategies are shifted towards increasing *c*_{Necrosis}.

### Data parameter space

Following the observation in [9], it was shown in-vitro that the proportion of neutrophils going into apoptosis or necrosis actually depends on the scale of insult. We utilize data on the peak values of necrotic population at time point 24 h. We also used an additional finding by Damas et al. where they show that 500 ng/ml LPS corresponds to the fatal concentration of LPS in humans, [31]. This corresponds to 100% necrosis, where the level of inflammation becomes uncontrollable due to the growing detrimental effects of local tissue damage generated from necrosis and the lack of anti-inflammatory functions from apoptosis. We are well aware that the assumptions we make using our model is limited by the data that we have at the moment. In particular, the power-law relation may turn out to actually be a different functional relationship when more data would be collected. However, the current work aims to provide baseline research to show what the concepts behind the model can do. Although the analytical derivations would change, we believe that the main methodology and qualitative conclusions would remain unchanged. Hopefully, our work motivates other researchers to take a deeper look into these concepts as well as invites stakeholders to generate additional data.

The data in [9] reports a 17% baseline level of necrosis even in absence of ITMs. This might be due to experimental setup and other laboratory procedures of the in-vitro experiment. However, it is reasonable to assume that no necrosis occurs in-vivo in absence of ITMs). Additionally, the data is not normalized since in the experiments not all neutrophils measured were in apoptosis or necrosis. For this reason, we correct the population of necrotic neutrophils in the data by normalizing it with respect to the total percentage of apoptotic and necrotic neutrophils used in the in-vitro study. This normalization is necessary since our model does not take into account the population of activated neutrophils that have not yet gone into apoptosis or necrosis. We show that, after such normalization, the percentage of necrotic neutrophils seem to exhibit a power-law behavior with respect to the initial concentration of ITMs that trigger an inflammatory response.

We model this power-law behavior in Eq. (10):

$$ p-\gamma =\beta \bullet {ITMs_{initial}}^k $$

(10)

where *ITMs*_{initial} corresponds to the initial concentration of ITMs, *k* = 0.0929, *β* = e^{−0.6} and a constant term *γ*, which we will explain shortly. In the log-log plot we observe a linear relation between percentage of necrosis with respect to the initial concentration of ITMs. This relation has an intercept on the y axis of −0.6 and a slope k equal to 0. 0929 (Fig. 4).

However, the data report 17% necrosis even without the presence of LPS. In order to remove this baseline amount of necrosis from the experimental dataset, we subtract a constant term *γ* = 0.17 in the power-law equation. We consider this baseline necrosis as originating from sources other than a systemic inflammation, which we do not take into account in our model.

In order to determine the parameter space of the EGT model, which reproduces this experimental relation, we compare Eq. (10) with Eq. (8). It is easy to see that *γ* and *β* are as follows:

$$ \boldsymbol{\gamma} =-\frac{\mathbf{\log}\left(\frac{{\boldsymbol{b}}_{\boldsymbol{Apoptosis}}+{\boldsymbol{c}}_{\boldsymbol{N}\boldsymbol{ecrosis}}}{\boldsymbol{\alpha} \left(\boldsymbol{m}-\boldsymbol{n}\right)}\right)+\boldsymbol{\alpha} \boldsymbol{n}{\boldsymbol{N}}_{\boldsymbol{T}}}{\boldsymbol{\alpha} \boldsymbol{m}{\boldsymbol{N}}_{\boldsymbol{T}}-\boldsymbol{\alpha} \boldsymbol{n}{\boldsymbol{N}}_{\boldsymbol{T}}} $$

(11)

$$ \boldsymbol{\beta} =\frac{\mathbf{1}\ }{{\boldsymbol{N}}_{\boldsymbol{T}}\left(\boldsymbol{m}-\boldsymbol{n}\right)} $$

(12)

We can reduce this 5-dimensional parameter space (*b*_{Apoptosis}*,c*_{Necrosis}*,m,n,α*) to 4 dimensions by extracting *n* from (12), which is the equation shown in (13).

$$ \boldsymbol{n}=\boldsymbol{m}-\frac{\mathbf{1}\ }{{\boldsymbol{N}}_{\boldsymbol{T}}\boldsymbol{\beta}} $$

(13)

Based on Eqs. (11) to (13), we are able to identify the parameter space for which the relation between necrosis and ITMs reported in the dataset is fulfilled. In order to visualize the parameter space that satisfies all our conditions we plot the solutions in Figs. 5, 6 and 7.

Nature, through a course of evolutionary time, has driven neutrophils to adopt these specific set of parameter values shown in Fig. 5. We show that the system evolved into stable percentages of neutrophil population that correspond to a combination of values of *α*, *b*_{Apoptosis}, *c*_{Necrosis}, where these three parameters have direct proportionality towards each other. That is, nature plays the evolutionary game and achieves stable percentages of neutrophil population by simultaneously opting for lower cost of going into necrosis while the global cost of remaining ITMs is also low. With higher cost of necrosis, higher values for the global cost of remaining ITMs is preferred. A higher cost for necrosis means further aggravating the current level of inflammation by inducing local tissue damage that increases the ITM levels in the system. The system, in response, imposes higher global cost of remaining ITMs in the system. This goes to show that nature might have made the necrosis death pathway purposely damaging to the organism so that the fitness will make sure that the events of necrosis are done only as few times as possible.

We visualize the same data parameter space in terms of parameters *m*, *b*_{Apoptosis}, *α* on the left panel (A) and *m*, *c*_{Necrosis}, *α* on the right panel (B) of Fig. 6. We show that when the threat of remaining ITMs is high, the evolutionary stable percentages of neutrophil population for high, or effective resolving capacity of neutrophils, correspond to a low range of values for benefits of apoptosis (A), and by Eq. (13), low range of values for costs of necrosis (B). This range of values for benefits of apoptosis and costs of necrosis become narrower and spans lower values when the threat of remaining ITMs is low. This is shown by the narrow wedge shape located at the bottom of the plot of panel A (shown in purple). This shows that nature could not have only made necrosis deliberately detrimental to the organism, it also could have made the anti-inflammatory benefits of apoptosis more effective especially when the presence of remaining ITMs in the system becomes increasingly threatening.

Furthermore, with stronger resolving power of ITMs, or the higher *m* is, the narrower the range of values for benefits of apoptosis and costs of necrosis, which spans lower values in the data parameter space. This goes to show how well the system compensates by choosing combinations of parameters that balance the intrinsic property of neutrophils to neutralize ITMs, the damaging effects of necrosis, and the anti-inflammatory benefit of apoptosis.

The preference for range of costs of necrosis and benefits of apoptosis at lower values with increasing ITM-resolving capacity *m* is also apparent in the data parameter space shown in Fig. 7. Nature could have made necrosis intentionally detrimental to the system, but when the intrinsic property of neutrophils to neutralize ITMs is favorable, costs for necrosis and the benefits of apoptosis adjust accordingly by having lower values. Nature could have chosen to maximize the organism’s fitness by limiting the number of neutrophils that go into apoptosis or necrosis when even a small amount of neutrophils can completely resolve the inflammation. By doing so, the organism has conserved a considerable amount of energy, which is favorable for its survival.

### Emergence of necrosis among neutrophils with limited interactions in cellular automata scheme

We use all the combinations of parameters derived in the mean-field scheme that corresponds to the experimental relation as input parameters for game played in the cellular automata algorithm. We do this under the assumption that these combinations of parameter values correspond to those that have been optimized during evolution by nature. The Global ITMs scheme assumes that an activated neutrophil bases its strategy on the strategy of its immediate neighbors and the total concentration of ITMs in the system.

Here we explore a lattice size of 50 × 50, so that the entire lattice is occupied by all 2500 neutrophils. In order to model systemic inflammation where ITMs are everywhere in the body, we 1) distribute the initial concentration of ITMs uniformly over the lattice and 2) assume periodic boundary conditions to assume a larger area of the tissue that corresponds to the entire body. Our results are summarized in Fig. 8.

The results using the Global ITMs scheme fits well with the data. On other hand, the percentage of necrosis in the Local ITMs scheme shows that regardless of initial ITM concentration, necrosis is always the dominant strategy and that this percentage of necrosis is maintained at about 100%.

The results of our simulations provide numerical evidence that indeed, neutrophils need sufficient information regarding the intensity of inflammation, which may biologically pertain to stimuli that originates directly or indirectly (via pro-inflammatory cytokines and other substances released by cells in the innate immune system) from the source of ITMs. On the other hand, relying solely on local stimuli that originate from surrounding cells would be an insufficient deciding measure for neutrophils to pick a death pathway.

Snapshots of what the Global ITMs lattice looks like for various initial concentrations of ITMs evolving through time are summarized in Fig. 9. A single time step in our simulation corresponds to a single iteration in the algorithm where an activated neutrophil chooses a strategy based on a computed fitness. It is apparent that in a system with lower concentration of ITMs, necrotic neutrophils are better off isolated. This implies that their payoffs are maximized if they are situated away from the other necrotic entities, most likely to avoid aggravating the cost of necrosis which we interpret biologically as local tissue damage. On the other hand, apoptotic neutrophils tend to form clusters and position themselves in between necrotic neutrophils to maximize their payoffs. However, the tendency of necrotic neutrophils to isolate themselves is not evident in a system with higher concentration of ITMs. At 500 ITMs, we observe that necrotic neutrophils tend to cluster together, which can be biologically interpreted as the process when necrotic neutrophils intentionally aggravate inflammation by inducing local tissue damage.

We emphasize that the simulations shown in Fig. 9 aim to show how certain strategies dominate and propagate in the system through the course of time. Showing these dynamics provide a snapshot of how the system behaves at a particular point in time based on rules we specified in the payoff matrix. With minimal insult (0 ITMs), the system opts to take apoptosis as the optimal choice of strategy. However, when faced with an intense magnitude of insult (500 ITMs), activated neutrophils prefer to go into necrosis *early on* in the iteration, assumedly as an attempt to purposely aggravate the inflammation, which, biologically speaking, enhances the innate immune system’s response by recruiting more immune cells into the tissue. Apoptosis comes on later in the iteration as an attempt to exercise their anti-inflammatory benefits. Indeed, the higher the magnitude of insult is, the higher the *overall* benefit necrosis contributes to the survival of the system, which is not only prominent at the end point of the iteration, but also is observed at the beginning of the iteration process. In fact, this corresponds well to what is observed in data [9]. The authors were able to show that through the course of 36 h, apoptosis slowly increases and dominates when the insult is minimal. However, with higher magnitude of ITMs, apoptosis remains minimal, and only shows a slight increase at the end of the experiment (see Fig. 2 in [9]).

Using cellular automata to model the choice of death pathway reveals behaviors that emerge from microscopic interactions. We investigate these microscopic behaviors by looking at the cumulative number of activated neutrophils that go into necrosis per unit timestep. Our results are summarized in Fig. 10.

The bold lines in the plots correspond to the average number of necrotic neutrophils obtained per timestep by setting 20 combinations of parameters which we have randomly chosen from the data parameter space. Here we show that the cumulative number of activated neutrophils going into necrosis stabilizes at points that follow a linear trend with respect to the iteration time step. Note that a single time step corresponds to a single iteration where an activated neutrophil chooses to go into apoptosis or necrosis. This is simply because a single iteration limits a single activated neutrophil to go into apoptosis or necrosis. In our simulation, the preferred strategy in the early part of the iteration is necrosis, which means that pro-inflammatory functions due to local tissue damage comes first, followed by apoptosis, which induces anti-inflammatory signals. Note that each iteration step here does not directly correspond to the biological time frame. That is, we only look at the end point where stable strategies are achieved.

### What kind of game are the neutrophils playing?

In this section we identify where the neutrophil game lies among existing games known in game theory. Based on the payoff values specified in Table 1, it is easy to see that what we have is a symmetric game, where the payoffs are only dependent on the strategies employed, but not on the players. More so, the payoff inequality of *D* > *C* > *B* > *A* perfectly describes a so-called deadlock game, where the strategy that is mutually beneficial is also the most dominant – apoptosis. However, due to the added effect of global cost to the average payoff as in (5), instead of having a pure strategy dominating as the Nash equilibrium (apoptosis, in this case), a mixed Nash equilibrium is observed as stable equilibrium.

It is easy to show that a mixed Nash equilibrium does not exist when the global cost term is not taken into account. The average payoff *F*(*q*, *p*) of an individual playing necrotic with probability “ ” in a population, which overall plays necrotic with probability “ *p* ” would then be given by:

$$ F\left(q,p\right)= qpA+q\left(1-p\right)B+\left(1-q\right) pC+\left(1-q\right)\left(1-p\right)D $$

(14)

We can calculate the optimal *q* by taking the derivative of the average fitness cost as shown below:

$$ \frac{dF\left(q,p\right)}{\delta q}= pA+\left(1-p\right)B- pC-\left(1-p\right)D=0 $$

(15)

Substituting the exact values from the payoff matrix in Table 1 leads us to the following inequality:

$$ -p{c}_{Necrosis}+\left(1-p\right)\left(-{c}_{Necrosis}+{b}_{Apoptosis}\right)-p{b}_{Apoptosis}-\left(1-p\right)2{b}_{Apoptosis}=0 $$

$$ -{c}_{Necrosis}-{b}_{Apoptosis}\ne 0 $$

(16)

The value *p* is undefined when searching for a mixed equilibrium, where the fitness of going into apoptosis and that of necrosis is equal. This shows that indeed, a mixed Nash equilibrium does not exist for the game of neutrophils when the global cost term is omitted from the average fitness equation. This also implies that the global cost is necessary to replicate the percentage of neutrophils observed in data. This exact mechanism also matches with the cellular automata scheme, where we show that information on the *global* concentration of remaining ITMs in the body via direct or indirect stimuli is necessary in order for the neutrophils to function as they do to resolve inflammation in the body. This global rule that appears to be the driving force for cooperation, which we refer to as *necrosis*, has been explored in a previous study [32]. Indeed, having global information of the system is sufficient to drive the percentages of necrotic and apoptotic population that are needed to resolve the inflammation in the body.