A mathematical model of combined CD8 T cell costimulation by 4-1BB (CD137) and OX40 (CD134) receptors

Combined agonist stimulation of the TNFR costimulatory receptors 4-1BB (CD137) and OX40(CD134) has been shown to generate supereffector CD8 T cells that clonally expand to greater levels, survive longer, and produce a greater quantity of cytokines compared to T cells stimulated with an agonist of either costimulatory receptor individually. In order to understand the mechanisms for this effect, we have created a mathematical model for the activation of the CD8 T cell intracellular signaling network by mono- or dual-costimulation. We show that supereffector status is generated via downstream interacting pathways that are activated upon engagement of both receptors, and in silico simulations of the model are supported by published experimental results. The model can thus be used to identify critical molecular targets of T cell dual-costimulation in the context of cancer immunotherapy.


S1 Stochastic multistate discrete logic modeling framework: background
The mathematical model, which takes on a stochastic multistate discrete logic framework, consists of 9 variables (Tr1, Tr2, Tr5, Nk, PKB, Bim, JNK, S, and C) and two external inputs (I and Ox) that can take on values in {0, 1, 2}, which correspond to low (0), medium (1), and high (2) activity, respectively. External input values are set to a constant throughout the duration of the simulation, and the variable values are initialized as described below. The system at time t can be represented with the state vector x(t) = {Tr1(t), Tr2(t), Tr5(t), Nk(t), PKB(t), Bim(t), JNK(t), S(t),C(t), I, Ox}, where x i ∈ {0, 1, 2} for all x i ∈ x. Then, the state vector updates in discrete time as follows, where Each f i is specified by the transition table for the variable x i (Figure 5), c is a continuity constraint which limits the maximum jump size for each x i to one level, and p i indicates the probability P that x i (t + 1) will update using g i , i.e.
We take p i = 0.5 for all i. We use an asynchronous update scheme, so that each variable is updated in a random order at each time step, therefore x(t * ) represents the state of the system x during the random update scheme, whereby a subset of x i ∈ x may have been updated. We note that the steady states for G(x(t)) are equivalent to those for G * (x(t)) = { f i } since steady states do not depend on the asynchronous, stochastic, or continuous update scheme in this framework [1].
We consider again the expository example presented in the Methods: the transition table for PKB is found in Figure  5(e), and we replicate it in Table S1. We take (Tr2, Ox, PKB)(t * ) = (0, 1, 2). Then, f PKB = 0. Due to the continuity constraint, g PKB = 1. Then, PKB(t + 1) = 1 if the system is updated (which it will be with probability p PKB = 0.5).

S2 Simulation of the model
We used Matlab R2017A to simulate the system. Initial conditions for the simulations can be found in Table S2. We set initial TRAF activity to low (at time t = 0, Tr1 = Tr2 = Tr5 = 0) since our focus is on 4-1BB and OX40-mediated increases in TRAF activity, and hence the downstream effectors are also set to either low (Nk = JNK = PKB = 0) or to moderate (Bim = 1) since Bim is repressed by TRAF1. We note that randomizing the initial conditions has no effect on the results (1000 simulations with random initial conditions were run, with a resulting zero variance for the steady state of all variables).

S2.1 Simulating knock-out and overexpression experiments
For ERK knock-out, we change the update rule for Bim inhibitors to only include PKB: if PKB is at its highest level (2), then we set Inh B = 1, else we set Inh B = 0. In this way, we remove the action of Tr1 towards Bim inhibition (since this inhibition is ERK dependent) and we do not allow maximal inhibition to occur even at high PKB levels since we assume that maximal Bim inhibition requires activation of ERK activity. For ERK overexpression, we set Inh B = 2.
In this way, we give a stronger weight to ERK-dependent mechanisms of Bim modulation over ERK-independent mechanisms. This is due to the presence of direct evidence that ERK is involved in 4-1BB (I)-mediated Bim inhibition [2]. We note that we do not model ERK directly as this would require us to distinguish between the strength of ERKdependent and independent Bim inhibition under normal physiological conditions of ERK. For PKB, JNK, and Nk over-expression and inhibition, we set the respective variables either to 2 (overexpression) or 0 (inhibition) throughout the duration of the simulation.