Divergent Resistance Mechanisms to Immunotherapy Explain Responses in Different Skin Cancers

Simple Summary Despite the successes of immune checkpoint therapy in treating metastatic skin cancers, most patients either fail to respond or become unresponsive to immunotherapy. We found that the cell–cell communication of two different immune cell types, macrophages and memory B cells, correlate very strongly with the response to immunotherapy. We built a mathematical model based on these results that predict different skin cancers would have different response rates and predict that a high ratio of memory B cells to macrophages is optimal for a response to immunotherapy. Abstract The advent of immune checkpoint therapy for metastatic skin cancer has greatly improved patient survival. However, most skin cancer patients are refractory to checkpoint therapy, and furthermore, the intra-immune cell signaling driving response to checkpoint therapy remains uncharacterized. When comparing the immune transcriptome in the tumor microenvironment of melanoma and basal cell carcinoma (BCC), we found that the presence of memory B cells and macrophages negatively correlate in both cancers when stratifying patients by their response, with memory B cells more present in responders. Moreover, inhibitory immune signaling mostly decreases in melanoma responders and increases in BCC responders. We further explored the relationships between macrophages, B cells and response to checkpoint therapy by developing a stochastic differential equation model which qualitatively agrees with the data analysis. Our model predicts BCC to be more refractory to checkpoint therapy than melanoma and predicts the best qualitative ratio of memory B cells and macrophages for successful treatment.

: Killing rate at which largest response is observed. For a fixed value, we determine the value at which the biggest response occurs. We see that melanoma experiences this maximal response at lower values of indicating that therapy is more likely to result in a response. By varying along the -axis, we see that this result is robust. To determine the size of the response we sought to maximize, we computed the rate of change of the largest, stable cancer burden as varied, normalizing for the size of the cancer.

Model reduction of the four-state model incorporating pro-inflammatory macrophages
Denote the population of pro-inflammatory macrophages as then we can derive a four-state model with

= (1 − ) −
(1) Here the additional terms account for the upregulation of anti-tumor immune cells and plastic conversion with anti-inflammatory macrophages by newly added pro-inflammatory macrophages.
We assume that the conversion within macrophages subpopulations is fast in timescale, hence we have the dynamical equilibrium = . Summing up the last two equations and utilizing this relationship, we then have the effective dynamics for anti-inflammatory macrophages which has exactly the same form with the M equations in three-state model except that the parameters are reformatted. In the effective B equations, the influence of additional term can also be equivalently reduced into cancer-regulated term , since they have the similar Hill-function term, and the dynamics of is positively regulated by . Therefore, the three-state model in main text can be viewed as the model reduction as four-state model, where the effect of pro-inflammatory macrophages has been accounted by the existing parameters.

Non-dimensionalization of the dynamical model and parameter selection
We non-dimensionalized our system to simplify the analysis but present our results in terms of these original equations. These are the non-dimensionalized equations (see Supplementary Table 3 for relationships of new variables to old variables):

Equilibria and their stability of deterministic model
In solving for equilibria, we are able to simplify the set of equations into a union of two solution sets: where is a degree 5 polynomial. Of the six roots, we select only the sensical ones, i.e. those where ( , , ) lies in the first octant of ℝ . We determine the stability of these fixed points using the eigenvalues of the Jacobian, looking for those with all eigenvalues having negative real part. In Figure 4B, we show the number of stable equilibria in the -plane. We choose to classify the stable fixed points based on the size of the cancer population. When = 0 is stable, this is complete elimination. When is nonzero but several orders of magnitude smaller than its carrying capacity, we describe this as a dormant state. When is at least 1% of the carrying capacity, we call this a high cancer state. When two such states are stable, we call the one with the larger cancer population very high.
Adding stochastic effect to the model We next turned to the reality of biological noise and considered how this could impact our model. Let be the state vector of our system at time and let ( ) be our time-independent ODE function. We take a generic, time-independent noise term, ( ) and generated the following SDE model: with being a standard Weiner process. For our noise term, , we assume there are both additive and multiplicative sources of noise for each population. These are given by independent Weiner processes and so is a 3 × 6 matrix. For determining the parameters of these functions, we chose them to be proportional to the parameters of our ODE system (see table below). For B cells and macrophages, these choices were tied to their source rates (additive coefficient) and proliferation rates (multiplicative coefficient). The energy landscapes shown in Figure 5A scaled these noise terms by 1/ . The matrix is (see Supplementary Table 4 for relationships of state variables and additive and multiplicative proportionality constants):