During each tumor time step Δt = 15min, for each tumor cell, a random tumor event is chosen based on the probabilities of proliferation and apoptosis. The probability of each event occurring during this time step follows an exponential distribution with a given rate of event. Tumor cells proliferate at a cell-dependent rate, which is the sum of a base rate α_T and an FGFR3-induced rate increase. This increase is directly proportional to the active FGFR3 dimer fractional occupancy ϕ_D , defined as the ratio of the concentration of active dimers on this tumor cell to the average concentration of total FGFR3 on tumor cells harboring the FGFR3 mutation. Proliferation of tumor cells is density-dependent, i.e. when the number of neighbors exceed a certain threshold, the tumor cell cannot proliferate. Moreover, tumor cells undergo apoptosis at a base rate δ_T . FGFR3 signaling decreases the rate of apoptosis and this decrease is dependent on ϕ_D of each tumor cell.

Tumor cells with low antigenicity have a fitness advantage compared to HA tumor cells in that LA tumor cells produce less ISF and are eliminated by CTLs at a slower rate. See Section 2.2 and Section S1 for further details.

A static vasculature model is included to model the influx of therapeutic agents and immune cells ( Section S1.3 ). Blood vessels are located on the border of the ABM lattice and lattice sites here are referred to as “perivascular”. Immune cells are recruited into the TME after each tumor update based on the size of the tumor at the start of the iteration. We assume that the rate of immune cells arriving in the TME is directly proportional to the tumor size. These new immune cells are placed randomly at empty perivascular lattice sites, from which they enter the TME.

To account for the faster timescale of immune cell migration, immune time steps are set to Δt_imm = 7.5min. At each time step, an immune event is randomly chosen from proliferation, apoptosis, movement, conjugation with a tumor cell, exhaustion and AICD, based on the probability of each event. The probability is calculated from the rate of each event in a similar way as tumor events.

Immune cells proliferate at a base rate of α_I unless the immune cell is either currently conjugated with a tumor cell or has already become exhausted. The proliferation rate is increased based on the local ISF concentration. As with tumor cells, immune cells must have sufficient space to proliferate. Immune cells undergo apoptosis at a base rate of δ_I at all times. If the CTL is engaged with a tumor cell when it is undergoing apoptosis, the CTL stops attacking the tumor cell. Non-exhausted immune cell move in the TME at a constant rate of movement, m. To allow for persistent movement in a single direction, and to improve simulation efficiency, immune cells move n _move steps at a time, with each step moving to a neighboring lattice site. The direction of movement is chosen randomly, but movement tends towards the direction with higher concentrations of ISF, accounting for the distance between lattice sites in the Moore neighborhood. Detailed calculations of movement gradient is found in Section S1.5.3 . Unengaged, active immune cells attempt to conjugate with a non-apoptotic, neighboring tumor cell at a constant rate, β. If the immune cell successfully engages the tumor cell, then the immune cell is labeled as engaged and starts to eliminate the tumor cell. We assume that immune cells employ the perforin/granzyme pathway to clear the HA cells, eliminating them in 30min (33). By contrast, immune cells use the Fas/FasL pathway to clear LA tumor cells, taking 2h to successfully induce apoptosis in the target cell. This difference in targeting mechanism follows from observations that in the absence of antigen, T cells preferentially employ FasL to target tumor cells (33).

Conjugation ends when either the tumor cell becomes apoptotic or the immune cell becomes exhausted. All immune cells in the model are assumed active upon reaching the TME and thus express PD-1 and are thus subject to PD-1 signaling, which can trigger exhaustion (34). The rate at which immune cells become exhausted is affected by the concentration of the PD-1-PD-L1 complex, following a Hill function. Exhausted immune cells wait to die and otherwise affect the system only by taking up space. Furthermore, immune cells can undergo AICD at a constant rate of d_a when they go long periods without conjugating with a tumor cell (35, 36).

To compute the amount of FGFR3 signaling and the effects of an FGFR3 inhibitor on tumor cells, we employ a global method developed in (37). Rather than using local concentrations of receptors, inhibitor, and complexes as state variables in an ordinary differential equation (ODE) for every tumor cell, we divide the TME into regions and update state variables averaged within these regions. To account for intra-region heterogeneity, we further divide each region into three subregions: non-mutantoccupied, mutant-occupied, and tumor-free. Section S1.6 contains full details of the system of ODEs describing FGFR3 dimerization, reactions between monomers, dimers and the FGFR3 inhibitor, diffusion of the inhibitor, as well as pharmacokinetics.

As discussed in Section 2.1, FGFR3 signaling alters tumor cell fate decisions by increasing the proliferation rate and decreasing the apoptosis rate. We also assume that FGFR3 signaling has downstream effects on the immune system. In accordance with observations that harboring an FGFR3 mutation correlates with lower CD8⁺ T cell infiltration (7), we assume that FGFR3 signaling decreases the immune recruitment rate by a factor dependent on the average ϕ_D value across all FGFR3 mutant tumor cells.

To determine the amount of PD-1 signaling on each immune cell, we make use of another implementation of a global method (37) similar to that used for FGFR3 inhibitor and a quasi-equilibrium assumption. We first solve reaction-diffusion equations for PD-1 inhibitor reacting with PD-1 on immune cells to obtain the average free PD-1 across all regions in the TME. This quantity is used as an initial condition for solving the PD-1-PD-L1 reaction to obtain the concentration of PD-1-PD-L1 complex, which determines the rate of exhaustion of immune cells as described in Section 2.2. Details of the equations are found in Section S1.8 .

We first analyze the effect of ICI on tumor growth and its efficacy’s dependence on the initial composition of the tumor. The initial FGFR3 mutant cell proportions are varied between 0% (wild type, WT), 50%, and 100% (mutant, Mut). The initial tumor antigenicity proportions are similarly varied between 100% low antigen (LA) cells, 100% high antigen (HA) cells, and a 50-50 split. At initialization, these features are assigned independently so that all included pairings are equally represented at the start of a simulation. We first observe that the 100% HA WT tumors ( Figure 2A , bottom-left) regress spontaneously even without treatment, indicating that at least one of these fitness advantages (loss of antigenicity or gain of FGFR3 mutation) must be acquired for progression. If only one is acquired, tumors grew, but ICI alone is successful ( Figure 2A , bottom row and left column). Importantly, this indicates that HA tumors retain sensitivity to ICI despite an FGFR3 mutation. We also note that in the LA WT case, ICI does eventually result in elimination, but only after the tumor nearly reaches carrying capacity. Finally, the remaining four panels represent tumors with a subpopulation of LA Mut cells, and none of these respond to ICI.

To understand the role of the immune response in effecting these outcomes, we looked at the CTL infiltrate throughout the TME at a time point prior to any of the observed peaks in tumor burden. We measure CTL infiltrate here as the percentage of all cells in the TME that are CTLs. In other words, before the tumor began to shrink. Thus, we selected Day 20 for the control arm and Day 12 for the ICI arm. In the control arm, only HA WT tumors regressed and contained more than 12% CTLs on Day 20 ( Figure 2B , top row). The dashed line indicates the 12% mark. In the ICI arm, however, the total CTL infiltrate was not an effective predictor of tumor response as no threshold could be drawn to divide responders and non-responders agnostic to initial tumor composition ( Figure 2B , bottom row). Response under ICI was more driven by antigen burden and absence of FGFR3 activation.

We next quantified the change in tumor composition under control and ICI. In both of these arms, the more fit cells (LA and Mut) gradually take over the tumor ( Figures 2C, D ). Under ICI, this shift accelerates so that the fitter, more immune-evasive cells compose more of the tumor at endpoint ( Figure 2D ). That is, the failure of ICI produces a tumor population with faster growth dynamics and more resistant to immune clearance.

To further understand the role of tumor composition on the efficacy of the immune response, we measured the spatial colocalization of CTLs within the tumor. Analogous to tissue sections, we first considered the density of CTLs within the middle z-slice of the tumor ( Supplementary Figure 1 ). We found that significant shifts in the CTL density occurred between the WT and the mixed mutant tumors ( Figure 3A ). Specifically, within a given antigenic status (columns of A), the comparison between blue (WT) and red (mixed mutant) always produces a significant difference in CTL density in both the control and ICI arms. Note we only show comparisons of single-step changes in the composition and within therapy arms. More specifically, we only show significant differences between neighboring colors of the same column or neighboring columns of the same color. We observe that this mutation-dependent pattern persists over time ( Figure 3B ) by computing the active CTL density in this convex hull throughout the simulation and grouping by therapy (rows) and antigen status (columns). The active CTL density in the absence of mutant tumor cells is consistently higher than that in the presence of mutant tumor cells, with the exception of the time period of tumor elimination observed in the rightmost column of Figure 3B .

To see if the immune activity was uniform throughout the tumor mass, we looked at the density of active and exhausted CTLs as a distance from the tumor center. By computing the probability density function (PDF) normalized by the volume of the spherical shell of each bin, we can identify the radii at which these immune cell phenotypes are enriched ( Figures 3C, D ). Note that as these are PDFs, their integral is 1, meaning these curves do not contain information about the total number of cells in each compartment. This allows a comparison between the relative enrichment on the same set of axes. The red curves in each panel show the tumor density, giving a baseline to compare against that is nearly uniform up to the leading edge of the tumor where this curve rapidly drops to 0. In the control case, both the active (blue) and exhausted (black) CTLs peak just inside the leading edge of the tumor and decrease towards the tumor center in most conditions. Under ICI, these peaks occur deeper in the tumor and the decrease in density towards the tumor center is less pronounced. This increased depth of penetration on ICI occurs despite the measurement occurring 8 days earlier than the control case, indicating that the CTLs are benefiting from ICI even far from the vasculature.

We next introduced a small molecule inhibitor of FGFR3 into the simulations to characterize potential synergies with ICI. To focus on the outcome of these simulations and make comparisons to mouse model experiments, we report the model metrics on Day 25, a typical endpoint for the mouse model experiments. Indeed, the in silico growth curves in Figure 2A show similar trends as our previously published mouse model experiments (18) ( Supplementary Figure 2 ). Using a Gaussian kernel to smooth the outcomes at Day 25, we see that anti-FGFR3 monotherapy does decrease Day 25 tumor burden for tumors with mutants present, as illustrated by the red peaks of the PDFs lying to the left of the black peaks in the middle and right column of Figure 4A . Nonetheless, the relative efficacy of anti-FGFR3 monotherapy compared with ICI depends on the antigenicity of the tumor, as seen by the different relative positions of blue and red peaks in the middle and right column of Figure 4A . Specifically, targeted therapy is most effective with LA tumors (top row) and least with HA tumors (bottom row). Measuring the efficacy of these therapies by their reduction in tumor cell count compared to control ( Figure 4B ), they exhibit synergistic effects, i.e., more than the sum of the individual effects, in tumors with LA mutants, i.e., the tumors that did not respond to ICI alone.

We then looked at the composition of the TME at the Day 25 endpoint. We first looked at the relative abundances of tumor subtypes and observed only modest shifts in composition across therapies conditioned on the initial composition ( Figure 4C ). Notably, there is a slight increase in the proportion of non-mutants (LA/blue and HA/yellow) under targeted therapy in the mixed mutant tumors (middle column). Regarding CTL infiltration into the tumor, the targeted therapy does increase the CTL colocalization with tumor cells, but only by a modest amount ( Figure 4D ). This helps explain the synergy between these two therapies: the anti-FGFR3 therapy neutralizes the proliferation and apoptosis advantages with little change in immune activity, while ICI increases the immune activity.

Having identified the synergy between these two drugs, we next test the sensitivity of this synergy to the FGFR3-mediated fitness advantages. We focus on the two cases in which we could achieve upwards of 30% reduction in tumor burden by Day 25: HA mixed mutants and HA mutants. We test the following four therapy schedules to compare against the control: FGFR3 monotherapy, ICI monotherapy, FGFR3 followed by ICI (FGFR3 1st), and ICI followed by FGFR3 (ICI 1st) ( Supplementary Figure 3 ). In the two combination therapies, the first therapeutic option is given in weeks 2-3 and the second is given in 3-4 ( Supplementary Figure 3 ). For each of these schedules, we test 50 parameter combinations of the FGFR3-related proliferation and apoptosis parameters. In Figure 5 , we display these on the x- and y-axes, respectively, by computing the proliferation rate and expected time to apoptosis assuming the FGFR3 dimerization reaction is at equilibrium without targeted therapy. For each of these 50 parameter combinations, we identify the minimal therapy that leads to the maximal response, which we define using a decision diagram ( Supplementary Figure 6 ). Briefly, we focus on a 30% reduction, i.e., some response, and a 90% reduction, i.e., a near complete response.

With a heterogeneous population with regards to the FGFR3 mutation, low proliferation rates of FGFR3 mutants results in a situation in which ICI monotherapy results in at least 90% reduction in tumor burden on Day 25 relative to control ( Figure 5A ). At higher proliferation rates, ICI monotherapy cannot produce even a 30% reduction ( Supplementary Figure 7A ). Instead, at a proliferation rate of 1.75d⁻¹, combination therapy sequenced so that ICI is given first can result in 30% or 90% tumor reduction when apoptosis occurs on the time scale of years or weeks, respectively. At proliferation rates above 1.75d⁻¹, these therapies are ineffective with one exception.

With HA mutants, the pattern is similar but with one notable difference. At lower proliferation rates, ICI monotherapy does produce a 30% reduction in tumor burden ( Supplementary Figure 7B ), but combination therapy is necessary to elicit a 90% reduction by Day 25 ( Figure 5B ). This is due to the longer timescale for ICI to reduce the tumor burden ( Supplementary Figure 2.1 ). Indeed, anti-FGFR3 monotherapy produces a stronger response initially due to its direct effect on tumor fitness and its faster pharmacokinetics ( Supplementary Figures 4 , 5 ). One consequence of this slower response to ICI is that the maximal tumor burden peaks at a high value. Even in the parameter regions in which ICI monotherapy results in 90% tumor reduction, this peak can be more than double the peak with ICI followed by targeted therapy ( Supplementary Figure 2.1 ).