Holistic View on the Structure of Immune Response: Petri Net Model

The simulation of immune response is a challenging task because quantitative data are scarce. Quantitative theoretical models either focus on specific cell–cell interactions or have to make assumptions about parameters. The broad variation of, e.g., the dimensions and abundance between lymph nodes as well as between individual patients hampers conclusive quantitative modeling. No theoretical model has been established representing a consensus on the set of major cellular processes involved in the immune response. In this paper, we apply the Petri net formalism to construct a semi-quantitative mathematical model of the lymph nodes. The model covers the major cellular processes of immune response and fulfills the formal requirements of Petri net models. The intention is to develop a model taking into account the viewpoints of experienced pathologists and computer scientists in the field of systems biology. In order to verify formal requirements, we discuss invariant properties and apply the asynchronous firing rule of a place/transition net. Twenty-five transition invariants cover the model, and each is assigned to a functional mode of the immune response. In simulations, the Petri net model describes the dynamic modes of the immune response, its adaption to antigens, and its loss of memory.


Introduction
In contrast to the human neural system, which is centrally coordinated by the brain, the immune system is decentralized. The lymph nodes are important compartments of the immune system. Every human organism has about 600 lymph nodes, which control and defend the different areas and organs, located at strategically important positions. Lymph nodes are perfused by lymph, which contains many substances and different forms of antigens. The lymph nodes show afferent and efferent lymphoid and blood vessels. One of the aims of the lymph nodes is to clean the body of internal and external antigens and apoptotic cells. The spatial differentiation of the lymph node and the concerted action of various cell types, e.g., "effector" T cells and "memory" B cells, are powerful strategies that enable the human body to survive the attacks of bacteria, viruses, and tumor cells.
The lymphatic system and immune response have been extensively studied by mathematical modeling. For an introduction and review, we refer to Margaris and Black [1] and Cappuccino et al. [2], respectively. Mathematical models involving the lymph nodes have covered various levels of abstraction, ranging from the molecular signaling pathways to whole-body models [3][4][5][6]. Appropriate modeling techniques include, among many others, PN provides a holistic view of the immune response of the lymphatic system. Quantitative models on specific processes may be superior in answering specific questions. For the development of a holistic quantitative model, however, the PN may be the best starting point. Figure 1. Sketch of a lymph node as a structured organ with compartments: blood vessels (red), subcapsular sinus (yellow), medulla (green), T zone (brown), and germinal center (light blue) with light zone (blue) and dark zone (dark blue). Lymph enters the subcapsular sinus via the topafferent vessel, flows around the lymphoid compartments, and may enter the T zone.
We tested approaches that had previously been conventionally applied only for molecular systems, such as metabolic systems and signal transduction systems. Here, we provide a biological interpretation of all invariant properties. We propose a workflow for systematically testing and verifying the network model during the construction phase. We formulate a set of structural requirements that a biologically meaningful PN model of cellular systems must fulfill. We perform an in silico knockout analysis of the model. Our simulations of the model were in accordance with the current understanding of the dynamics of an immune response, without using any kinetic data. The PN model enables further experiments to be run in order simulate hypotheses and generate new hypotheses.

Materials and Methods
A PN is a bipartite-directed weighted graph. The two types of nodes are called places and transitions; P and T denote finite sets of places and transitions, respectively. In our graphical representations (see Figure 2), places and transitions are represented by circles and squares, respectively. The set of arcs F ⊆ (P × T) ∪ (T × P) describes connections from places to transitions or vice versa. Each arc f ∈ F has a nonzero integer weight w( f ) ∈ N. Each place p ∈ P carries an integer number marking m(p) ∈ N 0 for the quantity of so-called tokens m 0 (p) ∈ N 0 , where m 0 is called the initial marking. In mathematical terms, a PN is defined as a quintuple N = (P, T, F, W, m 0 ) with: • a finite set of places P • a finite set of transitions T that are disjunctive from P, T ∩ P = ∅ • a set of arcs F ⊆ (P × T) ∪ (T × P) • the weights of arcs W : f ∈ F → w( f ) ∈ N • an initial marking m 0 : p ∈ P → m 0 (p) ∈ N 0 . Transitions denote cell interaction, activation, and migration. The lymph node (LN) is structured into the T zone (TZ, brown), germinal center (GC, blue), subcapsular sinus (SCS, yellow), medulla (ME, green), light zone (LZ, light blue), and dark zone (DZ, dark blue). Cells may enter and leave the lymph node via blood (BL, red) or tissue (TIS, not highlighted). For lists of places and transitions, refer to Tables A1 and A2 in Appendix A, respectively. For the reaction system, refer to Table A3 in Appendix A.
For each transition t, the set of arcs F defines a set of pre-places •t = {p | (p, t) ∈ F} and a set of post-places t• = {p | (t, p) ∈ F}. Border transitions are transitions without postplaces (output transitions) or pre-places (input transitions). Tokens are movable objects. The number of tokens on the places, i.e., the marking, can change according to a firing rule. Here, we apply the classical P/T (place/transition) firing rule. A transition is allowed to fire, i.e., is activated or has a concession, if and only if a sufficient number of tokens is located on all of its pre-places. The firing of transition t removes a number of tokens from each pre-place p ∈ •t according to the weight w( f ) of arc f = (p, t). Simultaneously, the firing of transition t adds tokens on each post-place p ∈ t• according to the weight w( f ) of arc f = (t, p).
A PN is called pure if no transition has a pre-place that is also a post-place. A PN is connected if an undirected path (ignoring the direction of arcs) exists between each pair of nodes. A PN is strongly connected if a directed path exists between each pair of nodes. A PN is called ordinary if all arcs have weights of one. A PN is conserved if the sum of tokens in all places does not change by the firing of transitions.
Each entry c ij of the incidence matrix C denotes the change of tokens to place i if a transition j fires. The number of occurrences of transition i in a firing sequence s of transitions is provided by the component x i = #i of a vector x : T → N 0 . In context-free grammars and Petri nets, the vector x is called the Parikh vector of sequence s. The minimal set of Parikh vectors x that satisfy the equation C × x = 0 is called TI. If all transitions are covered by TIs, the PN is called CTI. The minimal set of semi-positive integer vectors y that satisfy the equation y × C = 0 is called PI. If all places are covered by PIs, the PN is covered by place invariants (CPI). MIs are linear combinations of TIs that aim to be feasible for the empty initial marking. For a detailed definition and discussion, readers are referred to Amstein et al., 2017 [35].
We applied the software MonaLisa [42] to construct, visualize, and analyze our PN model. To determine various mathematical properties of the PN, we applied the INA tool [43,44]. For the in silico knockout, we applied the software isiKnock [33].

Results
We constructed a PN model of the lymph node. We focused on the interaction of cells in the compartments of the lymph node. Properties of the model such as PI, TI, and MI, were analyzed. In addition, we performed an in silico knockout analysis and simulations.

Construction
We constructed a PN model of cell migration, cell-cell interactions, cell activation, cell differentiation, and cell proliferation inside a lymph node. In the model, the lymph node (LN) is a structured organ with a subcapsular sinus (SCS), paracortex (TZ for T Zone), medulla (ME), and germinal center (GC) with light (LZ) and dark zone (DZ). For illustration and a list of abbreviations of compartments, refer to Figure 1 and Table 1, respectively. Inside the LN model, various cell types are responsible for the immune response: T cells, B cells, macrophages, follicular dendritic cells (DCs), and antigen-presenting cells (APCs). For a list of cells and their abbreviations, see Table 2. For an introduction to the architecture of LNs and cellular traffic in LNs, interested readers are referred to [45]. We constructed the PN model from scratch, starting with hand drawings by an expert pathologist in our group. The expert pathologist had decades of experience in diagnosis based on stained whole-slide images of LN and experimental research on the structure of LN, cell motility, cell-cell interaction, and cell properties; see [46][47][48][49][50][51][52][53][54][55][56][57][58][59] and citations therein. For a list of processes and their abbreviations, refer to Table 3. Our intention was to construct a basic model that considers only a minimal set of the major cellular players and compartments of the LN. The model was indented to mimic an adaptive immune response in accordance with the current understanding of the LN dynamicswithout using any kinetic data, as these are not available in most cases. For the analysis of the PN model of the LN, we constructed a PN without border transitions with an arc weights of one. Thus, the model represents a place-bordered and ordinary PN (PBOPN).  We applied the following two formal conditions for the construction of a biologically meaningful PN model of the lymph node: 1.
The PN is CTI. Each transition invariant has a specific biological meaning.

2.
Each place invariant of a PBOPN has a specific biological meaning. A detailed definition is provided below in Section 3.3.
Formal condition (1) is the completeness property of the model. Each transition has to contribute to the function of the model, i.e., has to be an active part of a functional biological module. TIs, well-known as elementary flux modes in systems biology, represent functional modules and pathways of biochemical systems [37,[60][61][62]. A transition that is not part of any TI is either unnecessary or indicates missing reactions [26,63,64]; hence, a biologically meaningful PN model should be CTI, expressing completeness.
Formal condition (2) guarantees a proper balance of the cells in the system. The total number of cells may change by cell death, cell replication, or the inflow and outflow of cells from and to the environment. We modeled the death, inflow, and outflow of cells by border transitions. Without border transition, processes such as death, inflow, and outflow would not be considered in the model. The PBOPN has no border transitions. To model cell replication, we introduced arcs with weights of two, e.g., reactions in the form of cell → 2 cells. The PBOPN has only arcs of weight one, and consequently lacks cell replication. In the PBOPN, the total number of cells is conserved. Special attention was required for the antibodies that are repeatedly produced by activated B cells. To keep the number of all species constant, we had to delete the production of antibodies. We required that the PBOPN without production of antibodies be conservative. This meant that each species had to be assigned to a PI describing its conservation. During the construction of the PN model, we found this check of the balance of cells and antibodies by the PBOPN to be worthwhile for the identification of biologically questionable or erroneous reactions. Figure 2 shows the full PN model of the lymph node. The PN has 49 places, 65 transitions, and 149 arcs. For a list of the places and transitions, along with their names and descriptions, refer to Tables A1 and A2, respectively, in Appendix A. Table A3 in Appendix A lists the complete reaction system. The file Data S1 in Supplementary Material represents the PN model in SBML (Systems Biology Markup Language) [65]. The background colors in Figure 2 indicate the compartments: T zone (brown), germinal center (blue), subcapsular sinus (yellow), medulla (green), light zone (light blue), dark zone (dark blue), blood (red), and tissue (white). For a list of compartments and their abbreviations, refer to Table 1. For the nomenclature of places and transitions, we applied an encoding based on the sub-strings in Table 2 and Table 3, respectively.

The Petri Net Model
In terms of PN formalism [43,44], the PN is pure, CTI, and connected. The PN has transitions without pre-places (input transitions) and transitions without post-places (output transitions). The PN has no places without pre-transition or post-transition. The PN is neither strongly connected, CPI, bounded, ordinary, or conserved.
The network describes the immune response of the lymph node to an influx of antigens. Antigen (AG) enters the system via the tissue (TIS) (see input transition IN_AG_TIS in the upper left part of Figure 2). Within the lymph, the antigen flows into the subcapsular sinus (SCS) of the lymph node (see the yellow-highlighted part of Figure 2). In the SCS, macrophages (M) recognize the exogenous invader, bind to the antigen, and translocate with the AG to the T zone (paracortex) (see the brown-highlighted part of Figure 2). In the T zone, antigen-presenting cells (APC) recognize the AG and activate T cells (T). Activated T cells bind to B cells and translocate to the light zone (LZ); see the light blue-highlighted part of Figure 2. In the light zone, the T cells are stripped off and activated B cells (B1), called centroblasts, translocate to the dark zone and proliferate [66,67] (see the dark bluehighlighted part of Figure 2). From the dark zone, the newly produced activated B cells (B2), known as plasma B ells, translocate via the germinal center (blue), T zone (brown), medulla (green), and blood (red) to the tissue. In the medulla (ME), plasma B cells can be stored as memory B cells (B3) and proliferate. In the tissue, the plasma B cells produce antibodies (AB). The produced AB specifically recognizes the AG and triggers the degradation of AG by macrophages in the tissue.

Place Invariants of a Place-Bordered Ordinary PN
To construct the PBOPN model, we set all arc weights to one and deleted all input and output transitions.  Table A7 in the Appendix A lists the ten PIs of the PN. These PIs describe the conservation of T cells (PI-0), antigen-presenting cells (PI-1), dendritic cells (PI-2), macrophages (PI-3), B cells (PI-4), activated cells (PI-5, PI-6), non-activated B cells, B cells with antigen (PI-7), and B cells with antigens PI-8 and PI-9.
The PN was not conservative, i.e., not covered by PI. Place AB_TIS was the only place not covered by a PI. Thus, the PN was not conservative. Tokens on place AB_TIS represented antigens that were newly produced by activated B cells in the tissue. Deletion of the arc to the place AB_TIS would have generated two additional PIs that represented the conservation of antibodies in the tissue.
Each PI of the PBOPN has a specific biological meaning. Each cell type could be assigned to at least one PI that described its conservation. No cells disappeared from the PN. Two pairs of PIs, (PI-5, PI-6) and (PI-8, PI-9), were identical except for a single place. Each pair described the conservation of a single species. The two redundancies of the PIs resulted from the functional coupling of species in certain reactions, e.g., the coupling between consumption of antigens and elimination of antibodies in the tissue; see transition SEP_M-AG-AB_TIS. Checking the cell balance by PI enabled control and discussion of the specific reactions for inflow and outflow, production, and elimination of species.

Invariants of the Full PN Model
In contrast to its PBOPN, the full PN had two input transitions and six output transitions. The two input transitions, IN_AG_TIS and IN_B_BL, described the infection of tissue with antigen and release of B cells from the bone marrow to the bloodstream, respectively. The six output transitions described the removal of species from the system, e.g., elimination of antigen, consumption of antibody, and death of B cells. The PN was not ordinary and had three arcs of weight two. One arc of weight two described the mitosis of activated B cells (B1) in the dark zone; see transition M_B1_DZ in the right bottom part of Figure 2.

Place Invariants
The full PN had six fewer PIs than its PBOPN. Antigens, cells that were activated by antigens, B cells, non-activated B cells, and B cells with antigens were no longer conserved. The four remaining PIs of the full PN describe the conservation of macrophages, T cells, dendritic cells, and antigen-presenting cells, respectively; see Table A4.

Transition Invariants
The full PN is CTI. Table A5 in the Appendix A lists the 25 TIs which cover the network. The entry in the "length" column provides the number of transitions of the TI. The majority of TIs (n = 13) are trivial, i.e., they contain only two transitions. The twelve non-trivial TIs contain between five and 28 transitions. As an example, Figure 3 highlights the largest TI, with 28 transitions. TI-8 represents the activation of a B cell by an AG, its migration from the blood to the para-cortex, and its subsequent passage via the light zone, dark zone, medulla, and blood to the tissue. Within the compartments, the B cell differentiates to a centroblast (place: B1_GC), a plasma B cell (place: B2_GC), and a memory B cell (place: B3_ME). In the pathway of TI-23, the differentiated B cell is neither replicated nor does it produce any AB. Without contributing to the adaptive immune response, the B cell is finally degraded as a plasma B cell in the tissue. TI-0, TI-10, and TI-23 are very similar pathways, each with small variations in the migration of the B cell. B cells produced by replication in either the dark zone (TI-16, TI-21, TI-22, TI-24) or the medulla (TI-14) may be degraded in the tissue. Additional degradation of activated B cells may take place in the dark zone. TI-5 corresponds to the steps from activation of B cells to degradation in the dark zone. TI-20 corresponds to the replication and degradation of a B cell in the dark zone. TI-9 corresponds to the production of AB in the tissue. The eight transitions of TI-9 represent the production of AB by plasma B cells in the tissue, recognition of the AG by the AB, and the destruction of the AG by macrophages; see Figure 4. to the T zone (place: B_TZ) and then, after activation, pass via the light zone, dark zone, germinal center, T zone, medulla, and blood to the tissue (places: B1_LZ, B1_DZ, B1_GC, B2_GC, B2_TZ, B2_ME, B3_ME, B2_BL, B2_TIS). B cells differentiate in several steps to become activated, and eventually become degraded in the tissue. For the abbreviations used to name the places, refer to Tables 1 and 2. For a list of TIs, refer to Table A5.
Note that each TI corresponds to a functional module of the immune response. However, the TIs do not offer a holistic view of the adaptive immune response. Only one TI, TI-9, represents the production of AB. The transitions of TI-9, however, include neither activation of B cells by invading AG nor any replication of activated B cells in the dark zone or the medulla. Essential couplings of the activation of B cells and the replication of activated B cells with the production of AB are not reflected by functional modules determined by TIs. A similar shortcoming of the concept of TI has been observed and discussed in relation to PN models of signal transduction systems. For PN models of signal transduction systems, the concept of MI has been shown to be advantageous [35].  Table A6  The majority of MIs in Table A6 in the Appendix A were impure. Only three MIs were pure: MI-2, MI-5, and MI-19. These three MIs describe the production and degradation of B cells in the blood in combination with the free movement of B cells to and from the T zone and germinal center. Neither an activation nor any interaction with other cells or antigens is involved in these pure MIs. This means that all activation processes and functional modules of the immune response are described by impure MIs. Each impure MI induces a subnet that contains at least one PI. The feasibility of impure MIs relies on the necessary condition of sufficient tokens on these PI; hence, none of the impure MIs can be feasible for a PN with empty initial marking.

Prediction of System Behavior via In Silico Knockout
Based on the MIs, we performed an in silico knockout analysis [36] and depicted the results in a knockout matrix. The knockout matrix in Figure 6 shows the influence of the invasion of antigen (transition: IN_AG_TIS), the production of B cells in the blood (transition: IN_B_BL), the replication of activated B cells in the dark zone (transition: M_B1_DZ), the replication of memory B cells in the medulla (transition: M_B3_ME), the production of antibody by activated B cells (transition: REL_AB_TIS) for the production of antibodies (place: AB_TIS), antigen (place: AG_TIS), and activated B cells with antigen in the light zone (place: B1_AG_LZ). The half-filled circle is colored red if the production of the corresponding species is no longer possible and green otherwise. Trivially, a stop of the invasion of antigen, i.e., a knockout of transition IN_AG_TIS, has the highest impact; antigen in the tissue (place: AG_TIS) and activated B cells with antigen (place: B1_AG_LZ) are no longer produced. Antibody (place: AB_TIS) production is, however, possible without antigens. A stop of the production of B cells, i.e., a knockout of transition IN_B_BL, forbids the production of newly activated B cells (place: B1_AG_LZ). A production of antibodies is, however, possible because of the replication of activated B cells and memory B cells in the dark zone and the medulla. A stop of the replication of activated B cells and memory B cells in the dark zone and the medulla, however, has no effect because each can compensate for the loss of the other. The only way to completely disable an immune response with antibodies would be a knockout of the production of antigens by activated B cells. The insensitivity of the antibody production to single knockouts of either antigen invasion, replication of activated B cells, or replication of memory B cells demonstrates the robustness of the immune response. Figure 6. Example knockout matrix. The name of a knocked out transition is indicated by the label of the row. The name of an affected place is indicated by the label of a column. Half-filled red circles denote that the knockout had an effect; green circles denote no effect of the knockout. The knockout matrix is computed by an in silico knockout based on Manatee invariants [36]. We chose five transitions and three places. The knockout of transition IN_AG_TIS (first row) stops the invasion of antigens. No antigen enters the tissue (red circle for AG_TIS) and no activated B cell enters the light zone (red circle for B1_AG_LZ), though the production of antibodies may continue to be active (green circle for AB_TIS). Overall, the majority of circles are green (no effect of a knockout). The knockout matrix indicates the robustness of the immune response to single knockouts; see text.

Prediction of System Behavior via Simulation
To predict the system's behavior in a deterministic way, we performed simulations of the token flow in the PN by applying an asynchronous firing rule. Figure 7A plots the time development of the tokens for the empty initial marking. Figure 7A plots the dynamics of antigens in the tissue (place: AG_TIS, green dots), B cells in the blood (place: B_BL, red dots), activated B cells in the dark zone (place: B1_DZ, yellow dots), memory B cells in the medulla (place: B3_ME, blue dots), and antibodies in the tissue, (place: AB_TIS, pink dots). Whereas the number of antigens increases linearly, the number of B cells in the blood fluctuates on a low level. Neither activated B cells in the dark zone nor memory B cells in the medulla are generated. There is no production of antibodies, and the system shows no immune response.
The lack of tokens on the four PIs hampers the functionality of the system; see Table A4 in the Appendix A. We started the simulation with 100 tokens on each of the places for antigen (AG_TIS), macrophages (M_TIS), antigen-presenting cells (APC_TZ), T cells (T_BL), and dendritic cells (DC_GC). Figure 7B plots the time development of the tokens for an initial marking. The simulation ran for 10,000 steps in an asynchronous simulation mode. The number of antigens (green dots) initially increases. After roughly 1500 steps, the systems adapts to the antigens. Activated B cells in the dark zone (yellow dots) and memory B cells in the medulla (blue dots) are produced, leading to an increase in the number of antibodies (pink dots). The number of antigens drops and fluctuates at a low level. The drop in antigen is a consequence of the increasing production of antibodies by activated B cells. The number of B cells in the blood (red dots) randomly fluctuates at a low level. The initial condition is an adapted system (10,000 steps) with no further inflow of antigens. After around 3500 steps, the number of memory B cells drops again to zero. The system loses its memory of the antigen. (D): The initial condition is an adapted system (20,000 steps) with no further inflow of antigens. A high number of memory B cells makes the system robust against memory loss. Even after 100,000 steps, the number of memory B cells remains nonzero. Eventually, the number of memory B cells randomly drops to zero and the system loses its memory to the antigen.
We saved the marking after 10,000 steps. We restarted the simulation with the saved marking, now without any inflow of antigens. Figure 7C shows the time development of the number of activated B cells in the dark zone (yellow dots) and memory B cells (blue dots). All activated B cells in the dark zone were lost after 1000 steps and were no longer able to reproduce. While the memory B cells persisted for a far higher number of simulation steps, they were eventually lost after around 4000 simulation steps. After 4000 simulation steps, the system had lost all activated B cells, and would have to adapt from scratch to a possible inflow of antigen. This simulation shows the memory loss of the immune response when not stimulated by the inflow of antigens.
We repeated the in silico experiments for a longer stimulation of the system by antigens. To this end, we saved the marking of the PN after 20,000 simulation steps with an inflow of antigens. Restarting the simulation with the saved marking and without an inflow of antigen yielded the time development in Figure 7D. Whereas the activated B cells in the dark zone (yellow dots) vanished after a few simulation steps, the number of memory B cells was nonzero for the entire simulation of 100,000 steps. Note that the number of memory B cells reaches low numbers several timesat around 20,000, 70,000, and 80,000 simulation steps. If the number of memory B cells had become zero at one of these simulation steps, the memory B cells would have been lost. The number of memory B cells is a result of replication and degradation. The random fluctuation induced by these two counteracting processes results in a certain non-zero probability that the number of memory B cells may reach zero. Without the simulation providing an inflow of antigens, the system eventually loses all memory B cells. After the loss of memory B cells, the system has to re-adapt to any new possible inflow of antigens. This simulation describes the finite time period of immunity to an antigen after an infection or vaccination.

Discussion and Conclusions
To predict the behavior of the immune system in the lymph node, we constructed a semi-quantitative PN model of the immune system. The components of the model were the major cellular players of the immune response. The advantage of the PN approach is an abstraction level that enables model construction, simulation, and verification even if knowledge of kinetic data, i.e., of numbers for reaction constants or the number of cells, are lacking. Our objective target was the development of a consistent model of the interplay of cells. To validate the consistency of the network, we applied a variety of theoretical methods and algorithms based on invariant properties.
Verification techniques based on invariant properties are well established for PN models of metabolic systems [68] and signal transduction [32,69]. Our model described motile cells that migrate between compartments of the lymph node, interact, replicate, and become activated. To the best of our knowledge, this is the first time that methods of theoretical verification based on invariant properties have been used and evaluated for a model of the concerted action of cells in the lymph node.
We proposed an iterative process of model construction and steps of verification to develop a consistent model. Our approach involved the incorporation of known biological processes and the analysis of invariant properties. We chose the abstraction level of the model with two objectives: first, to make the model as small as possible while including all of the processes that the expert medical pathologists in our group considered essential; second, from the theoretical side, we required that the network be CTI and that each transition invariant be assigned to a biologically meaningful process. Both of these conditions are well-known prerequisites of PN models of metabolic or signal transduction systems [26,32,63,64]. The CTI property indicates a complete and error-free model under the assumption of a system at steady state.
Additionally, we required that the token flow of the model fulfill biologically-inspired conservation properties. Note that a complex model structure may unintentionally describe a non-biological leakage or source of cells. Thus, each addition or removal of a species had to be explicitly modeled by either input/output reactions or replication cycles. To check the conservation properties of cells, we found it worthwhile to compute the PIs not of the network itself, but of the PBOPN. In the PBOPN, we deleted all input and output transitions and set the weights of all arcs to one. By setting all weights to one, we suppressed the replication of cells in our model. Note that other models may use other motifs to describe replication, and as such may require a different strategy to suppress replications. For both the ordinary and PBOPN cases, each place invariant had to correspond to the conservation of a specific cell type, e.g., B cells, macrophages, and T cells. The PBOPN had ten place invariants. In our model, we found one place, representing antibodies, that was not covered by a place invariant. The production of antibodies by activated plasma B cells was the biological justification for this non-covered place.
By applying the iterative process of model construction and verification of invariant properties, we developed our PN model of the lymph node. The final PN had 49 places, 65 transitions, and 149 arcs. The model was covered by 25 TIs and had four PIs. The four PIs represented the conservation of macrophages, T cells, dendritic cells, and antigen-presenting cells. The PBOPN had six additional place invariants representing B cells, activated B cells, and antigens. Note that the inflow of B cells and antigens was explicitly modeled by input transitions in the model. Activated B cells were produced by replication cycles.
The 25 TIs represented the diversity of functional modules. We assigned a specific biological function to each of the 25 TIs. However, we found essential couplings that were not reflected by TIs, e.g., no TI described coupling of the invasion of antigens to the activation of B cells, the production of antibodies, or the replication of activated B cells in the medulla. The shortcomings of TIs to describe complete pathways have been observed for PN models of signal transduction, motivating the concept of Manatee invariants [35]. Our PN model had 66 MIs. Among the 66 MIs, 14 MIs described the production of antibodies. The diversity of MIs describe the functional modules of the immune response in a holistic view, and include essential couplings, e.g., the coupling of antibody invasion to the production of antigen. Through an in silico knockout analysis, we found the production of antibodies to be quite robust against perturbations, such as the replication of activated B cells or the replication of memory B cells. A failure of persistent defense by the immune system would be lethal for an organism; hence, we expect any theoretical model of the immune system to be robust against perturbation.
We applied only known biological reactions. A list of the resources is provided in Table A2 in the Appendix A. To validate the combination of all reactions, we performed simulations. The simulation of our PN model predicted the dynamic behavior expected for an immune response. The systems adapted to infection by antigen and became able to control the number of antigens. The adapted model responded much faster to a second infection and was able to keep the number of antibodies at a low level. The effectiveness of the defense depended on the time period of the infection and on the time period between the first and second infections. After a long time period of the first infection, the defense mechanism became quite stable, though it eventually lost its memory to the antibody. After the loss of its memory, the model had to adapt anew in order to establish an effective immune response. The loss of memory was a random process that varied broadly from simulation to simulation. This behavior of the immune response is known from observations of the duration of the immune defense after infection or vaccination [70,71].
The ability of a PN model to mimic and predict the dynamic behavior of the immune response was surprising, as the model applied neither any quantitative rate constants nor any mass action kinetics. All reactions in the PN model had identical "time scales". No reaction was faster than another. We assumed that the majority of reactions of our model had similar rates and that the combinatorial diversity of pathways was the determining factor of the dynamic behavior.
Real-time imaging microscopy on human lymph nodes has recorded the velocities and 3D tracks of individual cells [72], while automated image analysis is able to measure the distribution of dimension of lymphatic compartments and the abundances of cells in the departments [57]. We plan to feed quantitative data into our PN model in order to develop a quantitative model. Eventually, we want to develop a model based on ordinary differential equations, stochastic simulation with the Gillespie algorithm, and an agent-based model.
In several aspects, this paper creates a unique basis for further investigations, especially clinical, pharmacological, and machine learning studies. The model presented here can be widened and adopted for the input of quantitative and qualitative human data. These strategies may enable us to understand deeper processes in human lymph nodes, as well as to predict complications and outcomes concerning CAR and NK cell therapies. In addition, despite the complexity of the system, the modulation of the cellular pathways in the model could elucidate aspects of cell interaction rules. Even without following primer causative suggestions, the results drawn from the input data may help in therapeutic decisions.   Table A1. The 49 places, their names, and a description of their biological meaning. The prefix of the name encodes the species, e.g., a cell, a set of interacting cells, or an antigen. The prefix is either a single abbreviation or a concatenation of abbreviations; for the abbreviations of types of species, refer to Table 2. The suffix of the name encodes the compartment in which the species is located; for the abbreviations of the compartments, refer to Table 1 M_AG_AB_TIS macrophages bind to the complex of antigen and antibody in tissue P-47 AG_dead_TIS antigen will be degraded in tissue P- 48 AB_used_TIS antibody, inactivated and marked for degradation (technical place) Table A2. List of 65 transitions, their names, and a description of their biological representatives. The prefix of the name encodes the reactions and processes. The abbreviations are explained in Table 3. In the middle are the interacting cells in the abbreviation from Table 2. The suffixes of the names encode the compartment in which the species is located; for the abbreviations of the compartments, refer to Table 1.

No. Name Description
T ENC_AG_M_with_APC_TZ Antigen-presenting cells contact the antigen and macrophages in T zone [94]  Macrophages contact the complex of antigen and antibody in tissue [92] T-64 SEP_M_AG_AB_TIS Antigens and antibodies will be marked by macrophages for degradation in the body (tissue) [92] Table A3 . Reaction system. The species "OExternal" denotes the interface to the environment. "OExternal" is ignored for preconditions and token balance.   Table A4. List of four place invariants (PIs). The columns "length", "places", and "conserved cell" provide the number of places of the PI, the names of all places of the PI, and the conservation of which cell (type) the PI originates from, respectively.  Activation of a B cell by an invading antigen. The activated B cell moves from blood to the T zone and via the light zone, dark zone, medulla, and blood to the tissue. The B cell differentiates in several steps and eventually becomes degraded in the tissue. Neither a cell replication nor the production of any antibody is involved. Activation of a B cell by an invading antigen. The activated B cell moves from blood to the T zone and via the light zone, dark zone, medulla, and blood to the tissue. The B cell differentiates in several steps and eventually becomes degraded in the tissue. Neither a cell replication nor the production of any antibody is involved. Activation of a B cell by an invading antigen. The activated B cell moves from blood to the T zone and via the light zone, dark zone, medulla, and blood to the tissue. The B cell differentiates in several steps and eventually becomes degraded in the tissue. Neither a cell replication nor the production of any antibody is involved. Activated B cells (B1) proliferate in DZ. B1 differentiates to plasma cells (B2) and moves via the light zone of the germinal center, T zone, medulla, and blood in the tissue. B2 is degraded in the tissue.

TI-22 12
ENC_B2_with_M_TIS DIF_B3_ME TRA_B3_ME_to_BL M_B1_DZ G1_B1_DZ SEP_M_B2_TIS OUT_B2_TIS TRA_B1_DZ_to_GC DIF_B2_GC TRA_B2_GC_to_TZ TRA_B2_TZ_to_ME TRA_B2_BL_to_TIS Activated B cells (B1) proliferate in DZ. B1 differentiates to plasma cells (B2) and moves via germinal center, and T zone in the medulla. In the medulla, B2 differentiates to memory B cells (B3). B2 moves then via blood in the tissue and will there be degraded. Activation of a B cell by an invading antigen. The activated B cell moves from blood to the T zone and via the light zone, dark zone, medulla, and blood to the tissue. The B cell differentiates in several steps and eventually becomes degraded in the tissue. Neither a cell replication nor the production of any antibody is involved.

TI-24 11
ENC_B2_with_M_TIS M_B1_DZ G1_B1_DZ SEP_M_B2_TIS OUT_B2_TIS TRA_B1_DZ_to_GC DIF_B2_GC TRA_B2_GC_to_TZ TRA_B2_TZ_to_ME TRA_B2_ME_to_BL TRA_B2_BL_to_TIS Activated B cells (B1) proliferate in DZ. B1 differentiates to plasma cells (B2) and moves via germinal center, T zone, medulla, and blood in the tissue. B2 is degraded in the tissue. Table A6. List of the 66 Manatee invariants (MIs). The columns "pure" and "antibody" indicate whether the MI is pure and whether the production of antibodies is part of the MI, respectively. The column "replication" indicates in which compartment(s) replication of activated B cells or memory B cells can occur; ME means the medulla and DZ means the dark zone. MIs that describe no replication of cells have an empty entry in the "replication" column.