4.2 Extension of the mathematical model
4.2.2 Law of large numbers/deterministic approximation and hybrid algorithm128
4 A stochastic model for melanoma T-cell therapy
of differentiated cells [NDif f/(NDif f +NDedi+NKO)] is weighed by a parameter α that determines the influence of this shielding effect. The factor (1−hNCyto)+ corresponds to cytokine-induced T-cell inhibition. As the number of inflammatory cytokines increases, the melanoma cells up-regulate PD-1 ligands (PD-L1) that inhibit the immune reaction. We comment on the necessity of adding the shielding effect and T-cell inhibition to the model in the next section, where we discuss the simulation results.
The ACTMETi therapy protocol is modelled through an addition of T-cells on the respective day of adoptive T-cell transfer and a temporary change of parameters during the five days of METi injections. These injections of METi are not modelled by an additional particle, which is reasonable since regular injections ensure a relatively stable level of METi, which is then degraded quickly after the injections stop. During the five days, melanoma cell reproduction, cytokine-induced dedifferentiation of WT melanoma cells, and T-cell inhibition are down-regulated, while reproduction and killing activity of T-cells are up-down-regulated, based on the experimental data in [87].
4.2.2 Law of large numbers/deterministic approximation and hybrid algorithm
4.2 Extension of the mathematical model
Assume that NK(0)/K → n0 ∈ Ω almost surely, as K → ∞. Then, for every T ≥ 0, (NK(t)/K)t∈[0,T] almost surely converges uniformly to(n(t))t∈[0,T], which is the unique solu-tion to the differential equasolu-tion
n(t) =˙ X
e∈E
veR˜e(n(t)) (4.6)
with initial condition n(0) =n0.
We want to apply this result to approximate the stochastic process that is described above with the deterministic solution to a differential equation. If, for some large but fixedK∈N, we can rewrite the generator ofN(t) with rates ˜Re as
Lφ(N) =X
e∈E
(φ(N+ve)−φ(N))Re(N) =X
e∈E
(φ(N+ve)−φ(N))KR˜e
N K
(4.7) and check the conditions of the theorem for these rates ˜Re, thenN(t)/K ≈n(t) with
n(t) =˙ X
e∈E
veR˜e(n(t)). (4.8)
To construct the new rates ˜Re, we have to rescale some of the parameters in the non-linear original rates Re with the fixed K. More precisely, we replace all competition parameters c∗,∗ with Kc∗,∗, the parameter for cytokine-induced dedifferentiation sCyDi,De with KsCyDi,De, the parameter for T-cell inhibition h with Kh, the parameter for T-cell proliferation bCD withKbCD, and the parameter for the killing of differentiated WT cellskDiwithKkDi. For example, if eis the event of cytokine-induced dedifferentiation,
Re(N) =sCyDi,DeNCytoNDif f =K(KsCyDi,De)NCyto K
NDif f
K , (4.9)
and hence ˜Re(n) = KsCyDi,DenCytonDif f. Apart from these changes, the new rates are the same as the original ones.
We check the conditions for these ˜Re on
Ω :=nn∈R{Dif f,Dedi,KO,Dead,CD8,Cyto}
+ :nDif f+nDedi+nKO > εo, (4.10)
for some small ε >0.
Condition (4.4) is automatically satisfied since E is finite, Ω0 is compact, and the only non-polynomial factornDif f/(nDif f+nDedi+nKO) is bounded by 1 for alln∈Ω.
To check condition (4.5), we note that, sinceEis finite, we can look at each eventeseparately.
We can also consider different factors individually, as long as they are bounded on compact subsets Ω0 ⊂Ω, since
ve( ˜R1e(n1) ˜R2e(n1)−R˜1e(n2) ˜R2e(n2))
≤ |ve|
R˜1e(n1) ·
R˜2e(n1)−R˜2e(n2) +
R˜2e(n2)) ·
R˜1e(n1)−R˜1e(n2)
. (4.11)
Since |(1−hn1,Cyto)+−(1−hn2,Cyto)+| ≤h|n1−n2|and 1/(nDif f+nDedi+nKO) is both bounded and Lipschitz-continuous away from nDif f +nDedi+nKO = 0, condition (4.5) is satisfied for our choice of Ω.
4 A stochastic model for melanoma T-cell therapy
Hybrid algorithm
We now want to apply this approximation result to speed up simulations for our stochastic model. In collaboration with the group of Prof. Martin Rumpf at the Institute for Numerical Simulations, in particular his student Kai Echelmeyer, we have developed a hybrid algorithm that combines a classical Gillespie algorithm with Runge-Kutta methods [66].
While a Gillespie algorithm produces an exact realisation of the stochastic process, it is computationally very heavy since every event is simulated separately. This particularly takes effect in large populations, where there are many reproduction and death events in a short period of time. In this case, the algorithm must perform many steps in order to simulate a time span of interest.
In contrast to this, the Runge-Kutta methods for the deterministic approximation of the stochastic process run much faster. However, they neglect any stochasticity, which is par-ticularly relevant in small populations. In those scenarios, random fluctuations could cause extinction of a subpopulation, while the deterministic solution will never quite hit zero. We want to keep these stochastic effects since for example the T-cells could die out and through this induce a relapse of the tumour.
In our algorithm we do not distinguish between small and large subpopulations, but rare and frequent events. We apply the deterministic approximation to the frequent events, i.e.
events with rates ˜Re above a certain threshold, which is not directly in line with the choice of Ω above. However, since the rates depend on the sizes of the different subpopulations, there is a direct correspondence and we can choose ε such that all population states that produce (partially) super-critical rates lie in Ω.
The stochastic-deterministic hybrid algorithm is implemented in C++. In the algorithm, the rates ˜R depend on the time as well as the state of the population. This is due to the fact that we allow for a change of parameters, e.g. during METi-injections. The threshold values lower(K) and upper(K), that determine whether an event is treated stochastically or deterministically, are chosen of orderK−1/2. This is because fluctuations of the rescaled stochastic process are of orderK−1/2 (compare for example the diffusion approximation in [171]) and we want to make sure that we allow for extinction of a subpopulation due to stochastic fluctuations. The parameters αij,βi, and γi determine the Runge-Kutta method and are in our case chosen as follows (summarised in the corresponding Butcher-tableau):
0
1 2
1 1 2
2 0 12
1 0 0 1
1 6
1 3
1 3
1 6
In the following we present the pseudo code for our stochastic-deterministic hybrid algorithm.
4.2 Extension of the mathematical model
More details can be found in [66].
Algorithm 1:Hybrid Algorithm Input: number of typesd, capacityK,
event set E with event vectorsve and rates ˜Re(t, n), stochastic eventsEs, deterministic eventsEd,
thresholds for updatelower(K)< upper(K),
maximal number of iterations kmax, end timeT, step size τ initial conditionn0
Output: Evolution of types {tk, n(tk)}k k= 0,t0 = 0,n(0) =n0
while k≤kmax andtk ≤T and Pdi=1ni(tk)>0 do fore∈ Es do
if R˜e(tk, n(tk))> upper(K)then Es=Es\e,Ed=Ed∪e
fore∈ Eddo
if R˜e(tk, n(tk))< lower(K) then Ed=Ed\e,Es=Es∪e
tk+1:=tk+τ,t←0,w←0 while t < τ do
R˜tot ←Pe∈E
s
R˜e(tk+t, n(tk) +w)) sample ttot ∼Exp(KR˜tot)
if t+ttot ≤τ then
sample e∗ ∈ Es proportional to{R˜e(tk+t, n(tk) +w)}e∈Es w←w+ vKe∗
t←t+ttot for1≤i≤4do
κi ←0 fore∈ Eddo
κi ←κi+ ˜Re
tk+γiτ, n(tk) +τPi−1j=1αijκj
n(tk+1) :=
K n(tk)+τP4 i=1βiκi
K
n(tk+1)←n(tk+1) +w k←k+ 1
4.2.3 Parameter choices
We have two sets of parameters corresponding to two sets of experiments. First, the exper-iments with pure wild type tumours from Glodde et al. [87], which are displayed in Figure 4.11A, and second, the experiments including PmelKO cells that are new to [88]. Between the sets of experiments, the melanoma cells showed slightly different behaviour (e.g. different speed of growth), however this can be achieved through variation of the parameters, leaving the systemic level (different events, structure of the rates) unchanged.
4 A stochastic model for melanoma T-cell therapy
In the following, we describe how the parameter values were either derived from the experi-ments to fit the model to the data or chosen to simulate scenarios beyond the experiexperi-ments.
The parameter choices for the two experimental setups and with or without the influence of METi injections are summarised in Table 4.1.
Initial cell numbers and scaling parameter K
We assume that the tumour takes the form of a 3-dimensional ball to relate the number of cells in our model to the tumour diameter measured in the experiments. A 5 mm tumour contains roughly 7·107 cells and since the measured tumour sizes vary from 1 mm≈5.6·105 cells to 10 mm ≈5.6·108 cells, we choose K = 107 as the typical size for the deterministic approximation. The number of initially injected melanoma cells is 2·105. Since likely not all of these cells survive and contribute to form the growing tumour, the initial condition NDif f(0) +NDedi(0) +NKO(0) is varied between 104 and 2·105.
Reproduction, death, and competition rates
The growth parameters of the different melanoma cells can be approximated fromin vivo ex-periments where tumour cells are injected into mice and then left to grow without treatment.
We set the death rates to 0.1 and assume the same rates for differentiated and dedifferenti-ated wild type tumour cells, since we cannot distinguish them in the experiments. Through a logistic fitting, we determine the birth and self-competition rates. The cross-competition between PmelKO and WT melanoma cells is only relevant to the new experiments. As mentioned in the previous section and seen in Figure 4.6, the growth rates of pure WT and PmelKOtumours are very similar and hence the competitive pressure that the WT melanoma cells impose on the PmelKO melanoma cells is chosen relatively large to fit the percentages in Figure 4.5B (results in Figure 4.16B).
Switch between differentiated and dedifferentiated melanoma cells
The rates for natural and cytokine-induced switching between differentiated and dediffer-entiated WT melanoma cells are chosen to fit experiments from Figure 4.2C, where WT melanoma cells are stimulatedin vitrowith TNF-αand/or METi for 72 hours, during which an equilibrium is attained. Without treatment, in a vehicle control experiment, a ratio of roughly 0.95:0.05 of differentiated to dedifferentiated melanoma cells is observed, which de-termines the ratio between the natural switch rates. The actual size of the rates is then chosen to make sure that an equilibrium is attained within 72 hours. This ratio shifts to approximately 0.65:0.35 under the influence of TNF-α, which indicates an additional dedif-ferentiation, induced by inflammatory cytokines. Under addition of METi alone, the ratio is about the same as in the control case, while a combination with TNF-α results in a ratio of roughly 0.85:0.15. This suggests that METi has no influence on the natural switch rates, while partially cancelling the cytokine-induced switch.
4.2 Extension of the mathematical model
T-cell activity and cytokines
The therapy parameters (reproduction, death, killing efficiency, and inhibition of T-cells, secretion and degradation of cytokines) are chosen to fit the experiments in Figure 4.11A and experiments with pure WT tumours corresponding to the protocol in Figure 4.1. During the course of therapy, 2·106 Pmel-1 T-cells are injected into the mice, which do not all infiltrate the tumour tissue. Simulations show that a variation of the number of T-cells has little to no effect on the evolution of the tumour since the amount of differentiated WT cells that can be killed is limited by the T-cell inhibition. Therefore, we fix the number of T-cells that are added on the respective day to 105 for all simulation runs. The variation between different mice is obtained by running simulations with different initial conditions (initial number of melanoma cells). The effect of METi on the proliferation, killing efficiency, and inhibition of T-cells is determined comparing experiments with therapy protocols including and excluding METi injections.
Dead melanoma cells
The parameter for the degradation of dead melanoma cells is chosen to fit the descent of the tumour size during the 5 days of METi injections. During this phase, the most dead melanoma cells are produced due to effective killing by T-cells, which contributes to the measured tumour diameter.
A
3 4 5
0 20 40 60 80 100
shielding parameterα
PmelKOindelfrequencyintumourcells(%) smallmedium
large
B
3 4 5
0 20 40 60 80 100
shielding parameterα
PmelKOindelfrequencyintumourcells(%) 19 mm10 mm
11 mm
Figure 4.8: (A) Percentage of PmelKO cells in simulations, taken at tumour size of 10 mm, varying parameter αand initial size of the tumour. Initial portion of PmelKO cells always at 17.1%. (B) Percentage of PmelKO cells in simulations, taken at varying tumour sizes and parametersα. Initial tumour of medium size with 17.1% PmelKO cells.
4 A stochastic model for melanoma T-cell therapy
Shielding effect
The parameter α, which scales the effect of differentiated cells being shielded from T-cell infiltration by other melanoma cells, is chosen to fit the percentages in Figure 4.5B (results in Figure 4.16B). The parameter has little influence in the simulations with pure WT tumours, which is why we determine it after fitting the other therapy parameters. An increase in α corresponds to a larger shielding effect. Since the experimental data in Figure 4.5B has a broad spectrum, we investigate several sources of variation to determine α. In Figure 4.8A, the composition of the tumour is analysed at a diameter of 10 mm, while the initial conditions are varied (keeping the initial percentage of PmelKOcells at 17,1%). Even though we compare initial sizes that lie below and above the critical threshold for therapy success, the percentages are very similar. Therefore, we fix the initial number of cells to a medium amount of 105 cells in Figure 4.8B and vary the time point at which the percentage is taken between the times when 9 and 11 mm diameter are reached. Since, in the experiments, the tumour is sequenced at the time when it has a size of at least 10 mm diameter, we setα to 4 to obtain an average of around 60% PmelKO cells between 10 mm and 11mm.
Table 4.1 summarises the parameter sets corresponding to the pure wild type experiments from [87] and the experiments with PmelKOvariants from [88]. Notice that these parameters correspond to the rescaled rates ˜Re and the deterministic differential equation. Therefore, some of them are rescaled byK (see the previous subsection for details).
[87], no METi [87], METi [88], no METi [88], METi
K 10e7 10e7 10e7 10e7
bDi=bDe 0.5 0.113 0.45 0.113
bKO - - 0.46 0.33
dDi =dDe=dKO 0.1 0.1 0.1 0.1
KcDi,Di =KcDi,De = KcDe,Di =KcDe,De
0.0053 0.0002 0.0047 0.0002
KcKO,KO - - 0.0048 0.003
KcDi,KO =cDe,KO - - 0.00002 0.00002
KcKO,Di =cKO,De - - 0.0117 0.0117
dDead 0.5 0.5 0.5 0.5
sDi,De 0.05 0.05 0.05 0.05
sDe,Di 0.95 0.95 0.95 0.95
KsCyDi,De 2 0.02 2 0.02
KbCD 60 70 70 100
dCD 0.4 0.4 0.4 0.4
KkDi 200 250 100 120
Kh 1.2 0.6 1.5 1
` 3 3 3 3
dCy 1 1 1 1
α 4 4 4 4
Table 4.1: Parameter choices for the different simulation scenarios.