Dynamic optimization elucidates higher-level pathogenicity strategies of Pseudomonas aeruginosa

Abstract

Multiple dangerous pathogens from the World Health Organization’s priority list possess a plethora of virulence components, including the ability to survive inside macrophages. Often, the pathogens rely on a multi-layered defence strategy in order to defend themselves against the immune system. Here, a minimal model is proposed to study such a strategy. By way of example, we consider the interaction between Pseudomonas aeruginosa and the human host, in which the host and the pathogen counter each other in a back-and-forth interaction. In particular, the pathogen attacks the host, macrophages of the host engulf the pathogen and reduce its access to glucose, the pathogen activates the glyoxylate shunt, which is started by the enzyme isocitrate lyase (Icl), the host inhibits it by itaconic acid, and the pathogen metabolizes itaconic acid using the enzyme succinyl-CoA:itaconate CoA transferase (Ict). The flux through the glyoxylate shunt allows the pathogen to avoid carbon loss and oxidative stress. These functions are of utmost importance inside a phagolysosome. Therefore, the pathogen needs to allocate its limited protein resource between the enzymes Icl and Ict in order to maximize the time integral of a flux through the enzyme Icl. We use both random search and dynamic optimization to identify the enzyme Ict as a cost-effective means of counter-counter-counter-defence and as a possible drug target during the early phase of infection.

dynamic optimization, nested defence, resource allocation, minimal model, Pseudomonas aeruginosa, counter-defence

Introduction

Pseudomonas aeruginosa is a widespread environmental microbe and a formidable human pathogen, which causes a range of infections that pose significant challenges to medical treatment (Stratton 1983, Mielko et al. 2019, Jurado-Martı́n et al. 2021). This Gram-negative, rod-shaped bacterium often affects individuals with compromised immune systems, such as those with cystic fibrosis, intensive care unit-acquired pneumonia, or chronic obstructive pulmonary disease (Fernández-Barat et al. 2017, Garcia-Nuñez et al. 2017, Camus et al. 2021). In the case of cystic fibrosis, the most widespread inherited lung disease worldwide, P. aeruginosa and Aspergillus fumigatus are the most common bacterial and fungal pathogens, respectively (Ferreira et al. 2015, Keown et al. 2020). Both of them cause other severe lung diseases (Emerson et al. 2002, Agarwal et al. 2020), and are involved in a complex interplay, which ranges from mutual suppression to cooperation (Briard et al. 2015, Ferreira et al. 2015, Reece et al. 2018). The metabolic byproducts of these pathogens can be harmful for the host (Dietrich et al. 2006, Gomez-Lopez et al. 2022).

One of the key features contributing to the resilience of P. aeruginosa is its intrinsic resistance to many antibiotics (Pang et al. 2019, Laborda et al. 2022). The bacterium possesses efflux pumps, multiple antibiotic-resistant genes including those of β-lactamases and the ability to produce biofilms, all of which contribute to its ability to resist conventional treatments (Sharma et al. 2014, Shapiro 2017, Lorusso et al. 2022), making it a persistent threat in diverse clinical settings (Mielko et al. 2019, ́Jurado-Martín et al. 2021). P. aeruginosa falls within the ‘critical’ category of the World Health Organization’s priority list of bacterial pathogens, demanding urgent research and development of new antibiotics (Tacconelli et al. 2018, Botelho et al. 2019).

An example of these adaptations is the ability of P. aeruginosa to survive and proliferate inside macrophages (Mittal et al. 2016). To that end, the bacterium utilizes the glyoxylate shunt, a pathway important for many pathogens (including Yersinia pestis, Mycobacterium tuberculosis, and the fungus Candida albicans) in order to survive inside eukaryotic host cells (McKinney et al. 2000, Anstrom and Remington 2006, Sasikaran et al. 2014, Riquelme and Prince 2020) and is, thus, considered a virulence factor. The glyoxylate shunt leads from isocitrate to malate (via glyoxylate) and succinate. The shunt allows bacteria to skip the steps of the tricarboxylic acid (TCA) cycle that produce CO₂ (Ahn et al. 2016). This simultaneously prevents carbon loss and reduces oxidative stress that is significant for a bacterium confronting with a macrophage (Riquelme and Prince 2020). That pathway also produces succinate, the substrate most readily consumed by P. aeruginosa during infection (Riquelme and Prince 2020), with lower protein costs. In order to produce succinate from isocitrate via the glyoxylate shunt, only one enzyme is needed rather than three in the classical TCA cycle. Moreover, the analysis of the elementary flux modes demonstrates that there are pathways leading to net production of succinate using the glyoxylate shunt (Schuster et al. 1999). At the same time, the shunt allows bacteria to utilize fatty acids as a carbon source (Lorenz and Fink 2002, Ahn et al. 2016).

In order to inactivate the glyoxylate shunt, macrophages produce itaconic acid, or itaconate. This metabolite is obtained by decarboxylation of cis-aconitate, one of the TCA cycle metabolites, and is a potent inhibitor of isocitrate lyase, the first enzyme of the glyoxylate shunt (Michelucci et al. 2013, Sasikaran et al. 2014). By this inhibitor, the ultimate level is not reached yet because to counter itaconate production, P. aeruginosa is capable of using an itaconate degradation pathway (Sasikaran et al. 2014).

Our understanding of the host–Pseudomonas interactions still lacks comprehensiveness; after 50 years of research, no anti-P. aeruginosa vaccines have been licensed yet (Sainz-Mejı́as et al. 2020). As will be shown in this paper, a theoretical framework called nested defence strategies, or defence and counter-defence (Schuster et al. 2019 Ewald et al. 2020) is useful for elucidating these interactions. Numerous cases of counter- and counter-counter-defence in nature were reported (Pumplin and Voinnet 2013, Hutin et al. 2015, Park et al. 2019, Hampton et al. 2020). In these terms, the chain of events in P. aeruginosa infections could be described as follows: the pathogen invades a host and invokes some harm (attack), the host’s macrophages react by engulfing the bacterial cells and putting them into phagolysosomes (defence). P. aeruginosa’s response is an adaptation to the new environment by utilizing the glyoxylate shunt (counter-defence). The macrophages, in turn, inhibit that pathway by producing itaconate (counter-counter-defence). In this paradigm, activation of the itaconate degradation pathway in P. aeruginosa is an example of counter-counter-counter- (or counter³-)defence. Such pathways have also been described for other pathogenic microorganisms, including species of the genera Micrococcus, Salmonella, Pseudomonas, and Yersinia (Sasikaran et al. 2014). Some enzymes of this pathway were also found in humans (Sasikaran et al. 2014) and M. tuberculosis (Wang et al. 2019).

In this paper, we use mathematical modelling to apply the nested defence strategy framework to the host–bacteria interaction. Mathematical modelling, coupled with clinical and experimental techniques, is an approach that was proven to be useful for understanding the processes of host–pathogen interactions and assigning optimal therapies (Waldherr et al. 2015, Dühring et al. 2017, Schleicher et al. 2017, Blickensdorf et al. 2019, Ewald et al. 2021). Modelling can help in devising wet laboratory experiments by highlighting the most promising paths and therefore saving researchers’ time and resources. Furthermore, some processes are hardly realizable experimentally, but can be readily simulated in silico (Chiacchio et al. 2014).

One of the techniques used in mathematical modelling of host–pathogen interactions is dynamic optimization, also known as optimal control (Lenhart and Workman 2007, Waldherr et al. 2015, Dühring et al. 2017). The underlying idea of dynamic optimization, an approach that initially arose from the field of engineering, is the necessity for living organisms not only to allocate their resources (metabolites, proteins, etc.) in the most efficient way, but also to be able to change this allocation in time to react to changing conditions (Klipp et al. 2002, Banga 2008, Bartl et al. 2010, 2013). In this framework, biological systems are usually described by ordinary differential equations, which include control (or decision) variables (Lenhart and Workman 2007, Banga 2008). These quantities, being also functions of time, are changed in order to optimize (minimize or maximize) an objective function.

Usually, models based on optimality principles in biology include side constraints. In several studies on resource allocation, a condition was used that constrains the sum of enzyme concentrations in a metabolic pathway to be constant (Klipp et al. 2002, Bartl et al. 2010, Ewald et al. 2024, Planque et al. 2024). The idea is that the cell can invest a certain protein capacity into the enzymes of the pathway. Here, we will include the constraint that the sum of the concentrations of two relevant enzymes either stays constant or changes in a predefined way because the cell is considered to increase its protein amount (capacity) allocated to the process in response to an external change, e.g. confrontation with the host.

Our main goal is to describe a host–pathogen interaction as an example of counter³-defence, using P. aeruginosa as a model organism and dynamic optimization as a suitable method of mathematical modelling. The resources of the pathogen are scarce while enclosed inside a phagolysosome, so the limited protein amount must be distributed between the pathogen’s enzymes in the most efficient way between counter- and counter³-defence. In the present paper, we set up a minimal model of the macrophage–P. aeruginosa interaction. Minimal models are useful when data about the system are limited and they can pinpoint areas for future investigation (Heinrich et al. 1977). Thereafter, we perform dynamic optimization in order to calculate the time courses of controls that maximize the objective function, namely, the flux through the pathogen’s glyoxylate shunt, and interpret the results.

Results

The model

Since our goal was to construct a minimal model, we only consider two species, three enzymes, and five metabolites (see Fig. 1). The following reactions are modelled: from the pathogen’s side, conversion of isocitrate to succinate and glyoxylate (ICL) and itaconyl-CoA formation (ICT); from the host’s side, only the reaction of itaconate production (CAD). The reactions are governed by the enzymes isocitrate lyase (Icl), succinyl-CoA:itaconate CoA transferase (Ict), and cis-aconitate decarboxylase (Cad), respectively. We assume Michaelis–Menten kinetics for the enzymes Ict and Cad. For the enzyme Icl, competitive inhibition by itaconate is used, based on observations (van Schaik et al. 2009). It is well known that the glyoxylate shunt operates in an irreversible way. Also, the other two reactions in the model can be considered irreversible (Michelucci et al. 2013, Sasikaran et al. 2014). For the bi-substrate reaction catalysed by Ict, the concentration of the co-substrate succinyl-CoA is considered constant and is integrated into the rate constant. The metabolites isocitrate and cis-aconitate are intermediates of the TCA cycle, which is thoroughly regulated in both pro- and eukaryotes (Millard et al. 2017, Arnold and Finley 2023). The enzyme Icl competes for its substrate isocitrate with isocitrate dehydrogenase (Icd), an enzyme of the TCA cycle. The turnover number of Icl was measured to be |$9.81$| s⁻¹ (Campos-Garcia et al. 2013), which is significantly lower than that of bacterial Icd (⁠|$85$| s⁻¹ ; Wang et al. 2019). In macrophages, the flux through the TCA cycle is diverted towards itaconate production; the new steady state is already reached before the simulation begins. Based on these facts, we can assume the concentrations of isocitrate (from the pathogen’s side) and cis-aconitate (from the macrophage’s side) to be constant, while the concentration of itaconate can vary. The fluxes through the system’s enzymes can be described by the rates of the corresponding reactions. Therefore, the model can be written as follows:

$$\begin{eqnarray} v\left( {Icl} \right) = \frac{{{k_{ca{t_{Icl}}}}\left[ {Icl} \right]\left[ {\textit{ICIT}} \right]}}{{{K_{{m_{Icl}}}}\left( {1 + \frac{{\left[ {\textit{ITAC}} \right]}}{{{K_i}}}} \right) + \left[ {\textit{ICIT}} \right]}} \end{eqnarray}$$

(1a)

$$\begin{eqnarray} v\left( {Cad} \right) = \frac{{{V_{ma{x_{Cad}}}}\left[ {Cad} \right]\left[ {CIS - ACO} \right]}}{{{K_{{m_{Cad}}}} + \left[ {CIS - ACO} \right]}} \end{eqnarray}$$

(1b)

$$\begin{eqnarray} v\left( {Ict} \right) = \frac{{{k_{ca{t_{Ict}}}}\left[ {Ict} \right]\left[ {\textit{ITAC}} \right]}}{{{K_{{m_{Ict}}}} + \left[ {\textit{ITAC}} \right]}} \end{eqnarray}$$

(1c)

$$\begin{eqnarray} \frac{{d\left[ {\textit{ITAC}} \right]}}{{dt}} = \frac{{{V_{ma{x_{Cad}}}}\left[ {Cad} \right]\left[ {CIS - ACO} \right]}}{{{K_{{m_{Cad}}}} + \left[ {CIS - ACO} \right]}} - \frac{{{k_{ca{t_{Ict}}}}\left[ {Ict} \right]\left[ {\textit{ITAC}} \right]}}{{{K_{{m_{Ict}}}} + \left[ {\textit{ITAC}} \right]}} \end{eqnarray}$$

(1d)

Figure 1.

Scheme of the proposed model. ITAC and CIS-ACO, host metabolites; ICIT, GLX and ITAC-COA, pathogen metabolites; Icl, Ict and Cad, enzymes; horizontal bar, regulatory influence, which connects ITAC with the glyoxylate shunt.

In Equations 1a–1d, the catalytic rate constants of Ict and Icl are indicated explicitly and multiplied by the respective enzyme concentrations because the latter are the control variables in the optimization. In contrast, [ITAC] is the state variable.

In the framework of this minimal model, the objective of the pathogen is understood as a maximization of the throughput of the glyoxylate shunt. Therefore, the objective function in our model is a time integral of the flux through Icl. Taking into account that the resources available for the pathogen inside a macrophage are scarce, we implement a constraint for the sum of the pathogen enzymes Icl and Ict. Similar constraints on the sum of enzyme concentrations have been used in optimality calculations earlier (Klipp et al. 2002, Bartl et al. 2010, Ewald et al. 2024). Therefore, the current optimization task can be formulated as follows:

$$\begin{eqnarray} \begin{array}{c} {\rm Maximize}\\ {\rm (Time\, integral\, of\, a\, flux\, through\, Icl)}\\ J = \smallint v\left( {Icl} \right)dt\end{array} \end{eqnarray}$$

(2a)

$$\begin{eqnarray} \begin{array}{c} {\rm subject\, to}\\ {\rm (Enzyme\, limitation)}\\ \left[ {Icl} \right] + \left[ {Ict} \right] = C \end{array} \end{eqnarray}$$

(2b)

The model parameters are collected or estimated from the literature; values and detailed descriptions are given in Supplementary Data S1.

In several simulations, we increase the capacity over time, representing the situation where the cell induces the defence after some stimulus. The change in capacity is a ‘decision’ by the cell that can be made very quickly. It is difficult to find an exact realistic value for that time span. We choose 80 s, which is in line with the finding that bacterial translation is very fast (Siller et al. 2010). This value is not critical, though, because the model is meant to explain the qualitative behaviour.

Model implementation

We model the situation when the pathogen is inside a phagolysosome and has already activated all the steps of its nested defence strategy. The k_catvalue of P. aeruginosa’s Ict was measured for a recombinant enzyme produced in Escherichia coli, and thus can differ from the one of an enzyme actually used by the pathogen. Therefore, we analyse three different scenarios according to different values of k_cat. These values may also arise by mutation or inhibition by some still unknown effector. We number these three scenarios by I, II, and III, with k_cat being equal to 100%, 10%, and 1%, respectively, of the value k_cat = 205.5 s⁻¹ (Sasikaran et al. 2014).

The steps we take to implement the model can be described as follows:

• Step 1. As a first approximation, we consider a simplified approach in which both concentrations of the enzymes Icl and Ict as well as the total protein capacity C remain constant throughout the time span of a simulation. The model under these assumptions is implemented in COPASI (‘COmplex PAthway SImulator’) (Hoops et al. 2006). A random search optimization with the objective function |$2{\rm{a}}$| is performed using in-built COPASI functions.

• Step 2. In order to estimate the benefits of using a counter³-defence mechanism for the pathogen, we perform the simulation in COPASI with the same parameter values as in Step 1, but with all the protein capacity allocated to the Icl enzyme. Thus, during this step we practically neglect the counter³-defence. Because of that, the value of the objective function will be the same for all three scenarios; this value is then compared with the results from Step 1.

• Step 3. Having a constant protein capacity is not biologically reasonable for a cell in many situations. Enzymes are costly to produce and to maintain; an evolutionary pressure exists not to produce enzymes that are not needed immediately. Neither the Icl nor the Ict enzyme is required by P. aeruginosa before confrontation with a macrophage. Therefore, we improve our model so that the protein capacity C is equal to zero at the beginning of the simulation, then increases linearly until it reaches the normalized value |$C = 1$|⁠, and remains constant afterwards. Then, we perform simulations of the time course of Icl and Ict using the enzyme allocation obtained in Step 1. To do so, the model is implemented in MATLAB. The resulting time course of enzyme concentrations for scenario III is demonstrated in Fig. 2.

• Step 4. We perform dynamic optimization of protein concentrations (considered as control variables), implementing the same time-dependent protein capacity as in Step 3. A special MATLAB toolbox based on a quasi-sequential approach (Bartl et al. 2011) is used. The basic concept of dynamic optimization applied to the host–pathogen interaction is illustrated in Fig. 3. Using our own work as an example, here the concentrations of metabolites such as itaconate are elements of the state variables vector |$x( t )$| and the concentrations of the enzymes such as Icl and Ict are elements of the control vector |$u( t )$|⁠. For a detailed explanation of the principles of dynamic optimization, see Ewald et al. (2020).

Figure 2.

The protein concentration time courses for scenario III (⁠|${k_{ca{t_{ICT}}}} = 2.055$| s⁻¹), calculated using dynamic protein capacity. Optimal protein distribution (see Table 1, column 3) was calculated by random search optimization for the case of constant protein capacity.

Figure 3.

The basic concept of dynamic optimization. The time courses are for illustratory purposes only and do not represent the current results. t, time; |$x( t )$|⁠, state variables (in our case, metabolite concentrations); |$u( t )$|⁠, controls (in our case, enzyme concentrations).

The introduction of a fixed time span in order to perform dynamic optimization leads to a boundary effect as can be observed in Fig. 4 for all of the scenarios: near the end of simulation, it becomes more beneficial to invest everything in Icl and neglect the itaconate metabolization by Ict. The negative effect of such a decision, namely a decreased level of Ict, does not have sufficient time to increase the concentration of itaconate and, thus, decrease the value of the objective function significantly. However, this behaviour is not biologically relevant. Therefore, we define a boundary area as the last |$5\% $| interval of the total time span and do not include the abrupt changes in the enzyme levels that had been obtained during this interval into the results presented in Table 1 (column 7). In order to be able to compare these results with the results of simulation without dynamic optimization, similar time limitations as for the dynamic optimization case are applied (Table 1, column 6). For additional information about simulation parameters and non-modified values of J, see Supplementary Data S2. The molecular masses of the enzymes Icl (UniProtKB entry Q9I0K4) and Ict (UniProtKB entry Q9I563) are |$58.887$| and |$43.079$| kDa. Nevertheless, we take the simple, unweighted sum of the two concentrations in the conservation quantity in Equation 2b in order to have a simplistic model. Moreover, properties other than the molecular mass also affect the enzyme costs, such as the amino acid composition.

Figure 4.

Enzyme concentration time courses obtained as a result of dynamic optimization for scenarios I–III (⁠|${k_{ca{t_{ICT}}}} = 205.5,\ 20.55,\ $| and |$2.055$| s⁻¹, respectively). The boundary effect is clearly seen for all the scenarios. Icl and Ict, concentrations of the corresponding enzymes; C, protein capacity.

Table 1.

Selected information about the simulation results for three different scenarios considered in this article.^a

Scenario	kcat, Icl, s−1	Protein allocated to Ict, (⁠|$C = 1$|⁠)	|$J\ ( {C = 1} )$|	|$J\ ( {[ {Ict} ] = 1} )$|	J (dynamic C)	J (dynamic optimization)
I	205.5	0.0448	347.2446	0.9378	249.7655	256.6999
II	20.55	0.1401	280.6844	0.9378	175.784	202.9297
III	2.055	0.4331	119.6270	0.9378	45.6222	84.0789

a

For the full table, see Supplementary Data S2. Column 4, results of random search optimization with the protein capacity being constant. Column 5, results of a simulation with protein capacity allocated solely to Icl. Column 6, results of the simulation using dynamically changing protein capacity and protein distribution found by random search optimization. Column 7, results of dynamic optimization with the protein capacity changing dynamically.

Ict synthesis is a cost-efficient mechanism of counter³-defence

We performed a random search optimization with constant protein capacity C for scenarios I, II, and III as described in Step 1 of the subsection ‘Model implementation’. The results are summarized in Table 1, column 3. Based on them we can conclude that, even in scenario III (the worst for the pathogen), no more than 50% of the protein capacity needs to be allocated to the Ict enzyme. The decrease in the objective function value, albeit significant (see Table 1, column 4), is not of the same order of magnitude as the decrease of k_cat. For all three scenarios, the allocation of all the protein C to the enzyme Icl (Step 2) leads to a very low objective function value, notably below unity (Table 1, column 5). Hence, it is two orders of magnitude lower than what can be achieved by using the Ict enzyme in addition (Table 1, column 4). These considerations allow us to identify Ict synthesis as a cost-efficient mechanism of the counter³-defence.

Dynamic regulation of resource allocation in the pathogen is only necessary for a slower Ict

We calculated the objective function values using the dynamic protein capacity as described in Step 3. These objective function values are collected in Table 1, column 6; an example of enzyme concentration time courses is demonstrated in Fig. 3. Our hypothesis was that, due to the trade-off between allocating protein capacity to either Icl or Ict, a dynamic regulation of such a resource allocation will be necessary. To test that hypothesis, we performed dynamic optimization of our model with dynamic protein capacity as described in Step 4 for scenarios I–III (Fig. 4); the results are summarized in Table 1, column 7. The comparison with the results obtained without dynamic optimization (Table 1, column 6) reveals that the difference between the objective function values for scenario I is ∼|$3{\rm{\% }}$|⁠. Therefore, in this scenario, there is no advantage for the pathogen to dynamically redistribute the fraction of protein capacity allocated to each enzyme by producing and maintaining a regulatory system capable of such a dynamic redistribution. It is also worth noting that the final enzyme distribution predicted by a dynamic optimization procedure is similar to the one predicted by random search optimization, the latter using constant protein capacity (compare Fig. 4 and Table 1, column 3). Therefore, constant protein capacity, albeit not biologically accurate, can in some cases be a decent approximation.

In contrast to scenario I, by comparing the results with and without dynamic optimization (Table 1, columns 6 and 7) for scenarios II and III, one can clearly see the possible benefits of regulating protein allocation. Using it in scenario III, it is possible to almost double the value of |${\rm{J}}$|. Therefore, it is very beneficial for a pathogen to have means of dynamic regulation for this system in case the k_cat value of Ict differs significantly from the one measured for the recombinant enzyme. However, the k_cat value of the recombinant Ict enzyme from Y. pestis was also measured to be relatively high (⁠|$375.98$| s⁻¹) (Sasikaran et al. 2014). That allows us to expect that the k_cat value of the wild-type Ict is closer to the one used in scenario I. If some mechanisms will be discovered in P. aeruginosa that perform such a dynamic protein reallocation, one can expect that macrophages possess some means to hinder a very productive Ict of the pathogen, and, therefore, not only a counter⁴-defence (Ict hindrance by host), but even a counter⁵-defence (Ict regulation by pathogen) may exist.

The obtained solution pattern is conservative in spite of parameter variations

Inspecting time profiles of the optimal solutions for scenarios I–III, one can observe one and the same strategy to be the most effective: First, all the protein capacity is allocated to Ict, then an abrupt switch occurs, and the new protein allocation remains relatively unperturbed until the end of the simulation. In order to test the stability of this strategy, we calculated log-normal distributions of the system’s kinetic parameters for scenario III (Ict’s |${k_{cat}} = 2.055$| s⁻¹), sampled 500 parameter sets, and performed dynamic optimization (Fig. 5). The same temporal behaviour pattern as in the previous simulations can be observed—not only for the best result, but also for all the conditions that lead to the |$50\% $| of the best results. Therefore, we can hypothesize that this pattern is a robust strategy that the pathogen should follow in case dynamic regulation of resource allocation is needed—e.g. when the Ict enzyme is hindered, or if its k_cat value is lower than estimated.

Figure 5.

The results of dynamic optimization for parameter sets composed from values obtained by log-normal distribution of kinetic parameters from scenario III. Icl and Ict, concentrations of the corresponding enzymes; |${\boldsymbol{C}}$|⁠, protein capacity; Icl with percentages in the subscript are the areas in which the best 50% and 10% of the solutions are located.

We have also performed sensitivity analyses based on three different approaches; see Supplementary Data S3 for details. The analysis indicates that there is a positive correlation between the value of the objective function, J, and the parameters |${K_{{i_{Icl}}}}$|⁠, |${k_{ca{t_{Icl}}}}$|⁠, |${K_{{m_{Cad}}}}$| and |${k_{ca{t_{Icl}}}}$|⁠. Between the J values and the remaining parameters, namely |${K_{{m_{Icl}}}}$|⁠, |${V_{ma{x_{Cad}}}}$| and |${K_{{m_{Ict}}}}$|⁠, there is a negative correlation. Moreover, the higher the parameter’s value, the smaller the impact of the parameter’s change (with a constant step width) on the non-normalized value of J. The only exception from this behaviour is the constant value of |$\partial J/\partial {k_{ca{t_{Icl}}}}$|⁠; this can be explained by the fact that, by definition, J is directly proportional to |${k_{ca{t_{Icl}}}}$|⁠. The absolute values of most response coefficients, |R|, are relatively small (less than unity) in the vicinity of the reference values, which indicates good robustness. The notable exceptions are the response coefficients for both |${k_{cat}}$| parameters in the model, due to the proportionality explained earlier.

Discussion

In this paper, we have treated an optimal resource allocation problem relevant for the pathogenic bacterium P. aeruginosa during a nested confrontation with its human host. To our knowledge, the presented minimal model is the first mathematical description of a counter³-defence in biology. In particular, the model describes the optimal choice of protein fractions invested into the enzyme isocitrate lyase (Icl), which initiates the glyoxylate shunt and in the enzyme succinyl-CoA:itaconate CoA transferase (Ict), which initiates itaconate metabolization. We have described this confrontation in terms of nested defence, formulated a minimal model, and studied it using the concept of dynamic optimization, also known as optimal control.

For analysing a counter³-defence in a qualitative way, a minimal model should be sufficient. Thus, we decided to consider the changes of only two enzyme concentrations (those of Ict and Icl). This does not exclude that other proteins are involved in the optimal allocation problem. In fact, by considering a gradual change in the capacity C, we did describe an addition or removal of protein from other enzymes.

We found that Ict production is a cost-effective means of counter³-defence, and that the given resource allocation problem does not require temporal regulation in order to be effective. The time course profiles obtained are robust towards parameter variations. We formulated a hypothesis that the putative temporal regulation of protein allocation in this system may indicate a presence of higher degrees of nested defence strategies.

The concept of dynamic optimization we applied to our model is currently used in various fields, including chemistry (Gerlinger et al. 2019), crime prevention (Ibrahim et al. 2023), epidemiology (Numfor et al. 2024), and systems biology (Ewald et al. 2017). Such a plethora of fields implies differences in scales of the models dynamic optimization is applied to: from protein complexes (Ewald et al. 2015, Soubrier et al. 2024) through metabolic pathways (Klipp et al. 2002, Bartl et al. 2010, Waldherr et al. 2015), like in the current work, to the population level (Numfor et al. 2024). The objective function can be formulated as a single function (Klipp et al. 2002, Ewald et al. 2017), or multiple functions (Waldherr et al. 2015), or even be separated between intra-organismal and population levels (Numfor et al. 2024).

In many cases, dynamic optimization problems can be solved by Pontryagin’s maximum principle in an analytical way, in particular if the system equations are bilinear in the sense that they are linear both with respect to state and control variables (Bartl et al. 2010, Ibrahim et al. 2023, Engida et al. 2024). Since our system is nonlinear in the state variable [|$ITAC$|], a solution by Pontryagin’s principle is very difficult if feasible at all. Therefore, we used a numerical solution procedure.

Interestingly, the optimal time profiles in form of rapid switches from one operating mode to another that we observed in this work were also reported by various authors that applied dynamic optimization to their models. For example, in one of the early works, Klipp et al. (2002) showed for a linear metabolic pathway that the optimal strategy will be to produce enzymes sequentially, with changes between expressed and inactive states as rapid as possible. Dühring et al. (2017) used dynamic optimization and game theory to study host–pathogen interactions, more specifically, non-lytic expulsion of pathogenic fungi by macrophages that engulfed them earlier, as a result of macrophages switching from an aggressive to a proliferative strategy. The approach of Waldherr et al. (2015) allowed them to describe known biphasic processes such as bacterial diauxic growth or usage of already excreted metabolites as a substrate. To further the understanding of fungal infection processes, Ewald et al. (2021) constructed a model of invasive aspergillosis. The authors illuminated the role of alveolar epithelial cells in the clearance of fungal spores and showed that the immune cell recruitment and depletion undergoes rapid changes throughout the infection course.

In our model, the proteins are assumed to be synthesized and degraded infinitely fast, which is, of course, not biologically realistic. However, such an idealized assumption was also previously used by Klipp et al. (2002) and Bartl et al. (2010). The problem has been discussed by Oyarzún et al. (2009), who modelled a metabolic pathway, first allowing infinitely fast switches. In a second simulation, they assumed enzyme production rates as controls instead, which led to finite slopes in the enzyme profiles. They showed that the time profiles of the solutions were similar in both cases. In the present model, we were mainly interested in finding a qualitative solution to the problem; therefore, we have adopted the idealized approach.

It is worth noting that the definition of a nested defence interaction and each of its steps is somewhat subjective. The first two steps of the P. aeruginosa–host interaction that is described in this work occur on a cellular level unlike the rest of them, which occur on a molecular level. If one considers the interaction of P. aeruginosa and macrophage (instead of the host), it will be an example of counter²-defence only, since the macrophage can be considered to be the attacking agent.

A further example where the classification is not clear-cut is the following: A bacterial pathogen attacks a host, the host produces a β-lactam antibiotic, the pathogen degrades this antibiotic by a β-lactamase, the host produces degradation inhibitors, and the pathogen tries to deactivate them as well. In this scheme involving a counter³-defence, clavulanic acid, a β-lactamase inhibitor that has, in addition, a modest antimicrobial activity (Filby et al. 2022), can play a role of defence (antibiotic) or counter-counter-defence (degradation inhibitor), or even both simultaneously. Moreover, although not occurring in nature, human-designed antibiotics such as sulbactam (also a β-lactamase inhibitor with additional antimicrobial activity; Penwell et al. 2015) or nanocarriers for clavulanic acid (Filby et al. 2022) are means to protect the host against pathogens. The addition of artificial nanocarriers may be an example of counter⁴-defence.

Theoretically, nested defence strategies can be built upon indefinitely by the competing organisms. In practice, the other examples of counter³-defence like the one described in this work were only scarcely reported (see earlier). One of the explanations could be that such a complicated system relies heavily on a presence of all its components, not only in the chosen organism, but also in its adversary. An inactivation of a defence level in an organism by a mutation or changing environment will turn the higher levels unnecessary (Schuster et al. 2019), and evolutionary pressure will occur in both organisms to stop maintaining the corresponding molecular mechanisms. On the other hand, when these mechanisms are gone, an evolutionary pressure may occur again, this time in order to encourage the invention of new nested defence layers. Therefore, the number of nested defence layers can be considered as a sort of trade-off between the necessity to outcompete the adversary and the effectiveness of cell resource management (Ewald et al. 2020). However, it should be noticed that the approaches of the host and the pathogen to this trade-off differ significantly: the host will most probably face more different pathogens during its lifetime than vice versa, and its genome is usually larger so that it can encode more layers of defence as well as utilize bet-hedging strategies more easily.

Genetic diversity allows P. aeruginosa to develop elaborated defence mechanisms depending on the environment, including the glyoxylate shunt and its initiation by the Icl enzyme, which is described in the current model. The glyoxylate shunt is important for a plethora of pathogenic species, including M. tuberculosis, in order to survive during the infection (Beste and McFadden 2010, Fahnoe et al. 2012, Sarkhel et al. 2022). An alternative mode of function of the glyoxylate shunt that is evolutionarily optimized to work under glucose limitation was predicted by theoretical methods (Liao et al. 1996, Schuster et al. 1999) and found later in experiment (Fischer and Sauer 2003, Beste and McFadden 2010). P. aeruginosa possesses the intricate regulatory system for activation/deactivation of the glyoxylate shunt that includes enzymes of both the systems that are present in M. tuberculosis and E. coli (Crousilles et al. 2018). Some microorganisms, like M. tuberculosis, have multiple variants of Icl; the necessity to counter all three of them is an insurmountable (up to date) obstacle on a road towards using this enzyme as a drug target (Antil and Gupta 2022).

Itaconic acid, or itaconate, is a small polyfunctional dicarboxylic acid. Being a potent inhibitor of Icl, itaconate even binds to it covalently in M. tuberculosis (Kwai et al. 2021). Itaconate is a core metabolite in our model, partaking in the nested defence interactions described by it. In general, itaconate is often mentioned as an important component of host innate immunity, having both anti-inflammatory and antimicrobial functions (Michelucci et al. 2013, Sasikaran et al. 2014, Ogger et al. 2020, Ku et al. 2022, Zeng et al. 2023). Recently, its usage for the treatment of inflammation-related heart diseases and non-alcoholic fatty liver disease was suggested (Ku et al. 2022, Weiss et al. 2023). Itaconate was shown to increase the permeability of P. aeruginosa biofilms for antibiotic tobramycin it forms complexes with (Ho et al. 2020). It is also used as a signalling molecule during P. aeruginosa pulmonary infection (Zeng et al. 2023). Some authors use itaconate as a structural part of nanocarriers for cancerostatics as well as antibiotics (also against P. aeruginosa) (Eid and Araby 2013, Kumar et al. 2023, Kumar and Kumar 2024). Interestingly, antimicrobial properties of the nanocarrier components, including itaconate monomers, were reported (Kumar and Kumar 2024).

Utilizing the itaconate degradation pathway that is initiated by the enzyme Ict, P. aeruginosa is able to turn the host’s itaconate production into an advantage for itself and use it as a carbon source (Riquelme et al. 2020). However, that requires a specific metabolic adaptation, to a degree that wild-type mice had significantly more severe infections than those unable to produce itaconate (Riquelme et al. 2020). In this paper, we focus on the early stages of infection, when such a metabolic adaptation has not yet taken place. Therefore, the flux through the glyoxylate shunt can be a suitable measure of fitness. It was also shown that, for the cases of idiopathic pulmonary fibrosis, the itaconate deficient mice, unlike their wild-type counterparts, develop persistent fibrosis (Ogger et al. 2020).

From our results, we can conclude that Ict as the counter³-defence is a cost-efficient way to secure the pathogen’s survival inside a phagolysosome and is, therefore, a suitable candidate for a drug target. Lowering the Ict enzyme effectiveness, as in scenarios II and III in comparison to scenario I, can lead to significant and partly opposing changes in a pathogen’s metabolism. Such an inhibition will render the pathogen unable to utilize the glyoxylate shunt due to the inhibition of Icl by itaconate. On the other hand, itaconyl-CoA, a product of Ict, inhibits bacterial methylmalonyl-CoA mutase (Mcm), an enzyme needed for propionyl-CoA-dependent growth (Peace et al. 2022). Ict inhibition could, therefore, increase the flux through Mcm, which produces succinyl-CoA, a metabolite of the TCA cycle. In view of side effects on host metabolism, it is worth noting that human Ict was shown to be homologous to the one from P. aeruginosa (Marcero et al. 2021). However, the corresponding gene was not shown to be expressed differently in macrophages during the interaction with P. aeruginosa, despite increased itaconate production (Yang et al. 2023). Therefore, it can be hypothesized that inhibition of bacterial Icl will influence host metabolism only weakly; this needs to be examined in detail upon drug development.

Many of P. aeruginosa’s virulence factors increase inflammatory activity such as influx of immune cells (examples being PQS, pyocyanin, PscI and PscF, ExoU) (Jurado-Martı́n et al. 2021). Since the bacterium can take advantage of immune cells and even facilitates their influx in order to do so, a favourable strategy for the host is to increase the bacterium's vulnerability towards immune cells. The aforementioned anti-inflammatory function of itaconate could theoretically help the host to counter some of the P. aeruginosa’s inflammation-related strategies.

Pathogenic microorganisms possess a vast array of virulence determinants to resist macrophages besides the counter-defence strategy analysed here. P. aeruginosa can inhibit phagocytosis using the distal parts of lipopolysaccharides (OPS) protruding from the outer membrane (Huszczynski et al. 2019) or use the porine OprF to avoid destruction in phagolysosomes, albeit by an unknown mechanism (Moussouni et al. 2021). It has been suggested for the fungal pathogen Candida auris that the host’s immune system is overcome by outcompeting macrophages for glucose and eventually starving them (Weerasinghe et al. 2023). C. albicans can utilize the glyoxylate shunt in order to facilitate gluconeogenic growth and eventually escape macrophages (Lorenz et al. 2004). However, it is unclear whether glucose starvation of macrophages is a relevant scenario in vivo and whether these mechanisms are also relevant for bacteria. The aforementioned virulence determinants, along many others, can also theoretically be countered by the host, thus giving rise to a plethora of possible nested defence strategies worth future studies.

By our analysis, we have identified the enzyme Ict as a possible drug target in order to combat P. aeruginosa infections. The present model can also be adapted to other pathogenic species that rely on itaconate metabolization during a confrontation with a host, e.g. the plague bacterium Y. pestis. Taking into account the cooperative nature of some interactions between P. aeruginosa and A. fumigatus during their co-infection, our work can indirectly help to treat A. fumigatus infections as well. However, experimental studies are needed in order to estimate the clinical benefits and relevance of the suggested drug target.

Acknowledgements

The authors thank Jan Ewald for helpful discussions and support with the optimization toolbox.

Author contributions

W.D. and S.S. conceptualized the study and established the mathematical model. W.D. collected the data, performed numerical simulations, and produced the figures. S.S. coordinated and supervised the study. Both authors wrote and reviewed the manuscript.

Conflict of interest

None declared.

Funding

W.D. and S.S. gratefully acknowledge financial support by the German Research Foundation (DFG) through the TRR 124302 FungiNet, ‘Pathogenic fungi and their human host: Networks of Interaction’, project number 210879364. We acknowledge support by the German Research Foundation (project number 512648189) and the Open Access Publication Fund of the Thüringer Universitäts- und Landesbibliothek Jena.

Data availability

The code used for modelling and the instructions on how to reproduce the results can be found at https://git.rz.uni-jena.de/qe45cow/p-aeruginosa-counter3-defence.