Next Article in Journal
Asymmetric Information and Credit Rationing in a Model of Search
Next Article in Special Issue
High Cost of Survival Promotes the Evolution of Cooperation
Previous Article in Journal
Two-Valued Strongly Group Strategy-Proof Social Choice Functions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Computing Stackelberg Equilibrium for Cancer Treatment

Ganzfried Research, Miami Beach, FL 33139, USA
Games 2024, 15(6), 45; https://doi.org/10.3390/g15060045
Submission received: 8 October 2024 / Revised: 29 November 2024 / Accepted: 18 December 2024 / Published: 23 December 2024
(This article belongs to the Special Issue Evolution of Cooperation and Evolutionary Game Theory)

Abstract

:
Recent work by Kleshnina et al. has presented a Stackelberg evolutionary game model in which the Stackelberg equilibrium strategy for the leading player corresponds to the optimal cancer treatment. We present an approach that is able to quickly and accurately solve the model presented in that work.

1. Introduction and Problem Formulation

In this section we review the Stackelberg evolutionary game dynamic model of cancer evolution previously studied [1]. There are two players: a follower and a leader. The leader is a physician who selects amounts of two different drugs to use for therapy, m 1 and m 2 . The follower is a cancer population consisting of three cell types: 0-type denotes a cell that does not develop resistance to either drug, 1-type cells are resistant to just drug 1, and 2-type cells are resistant to just drug 2. The follower selects a population size for each type, denoted x 0 , x 1 , x 2 , as well as a trait for each type, denoted u 0 , u 1 , u 2 . It is assumed that each of these variables are nonnegative, with the u i [ 0 , 1 ] . It is also assumed that they are all implicitly functions of time t . Note that u 0 does not appear anywhere in the analysis so can be ignored.
For each cell type i, there is a fitness function G i ( u i , m , x ) that the follower is trying to maximize. We assume that the dynamics of the population x are governed by
x ˙ = G ( t ) x .
In order to ensure that we are in equilibrium of the ecological dynamics we must have that x i ˙ = 0 for i = 0 , 1 , 2 . Thus, the follower is selecting u 1 that maximizes G 1 , u 2 that maximizes G 2 , and x 0 , x 1 , x 2 that ensure equilibrium of the ecological dynamics. The leader, knowing that the follower will subsequently select their actions in this way, selects m 1 , m 2 to maximize a quality of life function Q ( m , u , x ) .
Thus, we can formulate the problem of determining the optimal strategies for both players as follows:
max m * , u * , x * Q ( m * , u * , x * ) s . t . x ˙ i * = 0   for   i = 0 , 1 , 2 u i * arg max u i G i ( u i , m * , x * )   for   i = 1 , 2 m * 0 , x * 0 , 0 u * 1
This general model is instantiated by the following functional forms for the fitness functions G 0 , G 1 , G 2 and quality of life function Q. Note that the model presentation is slightly different between the paper [1] and the implementation in the code repository [2]. We will be using the model presented in the code.
G 0 = r max 1 α 00 x 0 + α 01 x 1 + α 02 x 2 K d m 1 k 1 m 2 k 2
G 1 = r max e g 1 u 1 1 α 10 x 0 + α 11 x 1 + α 12 x 2 K d m 1 b 1 u 1 + k 1 m 2 k 2
G 2 = r max e g 2 u 2 1 α 20 x 0 + α 21 x 1 + α 22 x 2 K d m 1 k 1 m 2 b 2 u 2 + k 2
Q = Q max c x 0 + x 1 + x 2 K 2 w 1 m 1 2 w 2 m 2 2 r 1 u 1 2 r 2 u 2 2
The model has several parameters, whose interpretations are summarized in Table 1. Note that in the code additional parameters a 0 , a 1 , a 2 , a 3 are defined, with
α 00 α 01 α 02 α 10 α 11 α 12 α 20 α 21 α 22 = a 0 a 1 a 1 a 2 a 0 a 3 a 2 a 3 a 0

2. Prior Approach

In this section we present the approach described in Github repository created by the authors [2]. They first solve for expressions for x i in terms of u i , m i such that the condition for equilibrium of the ecological dynamics is satisfied. This involves solving a system of three equations G i x i = 0 with unknowns x 1 , x 2 , x 3 . They calculate the following analytical solution:
x 0 * = K ( X 01 + X 02 + X 03 ) r max ( a 0 a 3 ) ( a 0 2 2 a 1 a 2 + a 0 a 3 ) ,
where X 01 , X 02 , X 03 denote the following quantities:
X 01 = a 0 2 a 3 2 d m 1 k 1 m 2 k 2 + r max
X 02 = a 1 ( a 3 a 0 ) e g 1 u 1 d m 2 k 2 + e g 1 u 1 r max m 1 k 1 + b 1 u 1
X 03 = a 1 ( a 3 a 0 ) e g 2 u 2 d m 1 k 1 + e g 2 u 2 r max m 2 k 2 + b 2 u 2
x 1 * = K ( X 11 + X 12 + X 13 ) r max ( a 0 a 3 ) ( a 0 2 2 a 1 a 2 + a 0 a 3 )
X 11 = a 2 ( a 3 a 0 ) d m 1 k 1 m 2 k 2 + r max
X 12 = ( a 0 2 a 1 a 2 ) e g 1 u 1 d m 2 k 2 + e g 1 u 1 r max m 1 k 1 + b 1 u 1
X 13 = ( a 1 a 2 a 0 a 3 ) e g 2 u 2 d m 1 k 1 + e g 2 u 2 r max m 2 k 2 + b 2 u 2
x 2 * = K ( X 21 + X 22 + X 23 ) r max ( a 0 a 3 ) ( a 0 2 2 a 1 a 2 + a 0 a 3 )
X 21 = a 2 ( a 3 a 0 ) d m 1 k 1 m 2 k 2 + r max
X 22 = ( a 1 a 2 a 0 a 3 ) e g 1 u 1 d m 2 k 2 + e g 1 u 1 r max m 1 k 1 + b 1 u 1
X 23 = ( a 0 2 a 1 a 2 ) e g 2 u 2 d m 1 k 1 + e g 2 u 2 r max m 2 k 2 + b 2 u 2
While it is not given in the model formulation, the code assumes that m 1 and m 2 fall in the interval [ m ̲ , m ¯ ] , with m ̲ = 0 , m ¯ = 1 . They discretize the space so that m 1 and m 2 are multiples of h = 0.1 , and iterate over all combinations, of which there are 1 h 2 . For each combination of values they calculate the optimal values of u 1 , u 2 as follows. They assume that the optimal value of u 1 is where G 1 u 1 = 0 , and the optimal value of u 2 is where G 2 u 2 = 0 . They calculate the following expressions for the partial derivatives:
G 1 u 1 = g 1 r max e g 1 u 1 1 α 10 x 0 + α 11 x 1 + α 12 x 2 K + m 1 b 1 ( k 1 + b 1 u 1 ) 2
G 2 u 2 = g 2 r max e g 2 u 2 1 α 10 x 0 + α 11 x 1 + α 12 x 2 K + m 2 b 2 ( k 2 + b 2 u 2 ) 2
To find values of u 1 , u 2 for which these derivatives equal zero, they perform an optimization to minimize the quantity | G 1 u 1 |   +   | G 2 u 2 | . They perform this optimization in Matlab using the fminsearchbnd function, using bounds 0 u 1 1 , 0 u 2 1 . Note that the above expressions for x i * which are in terms of m i and u i are substituted into the expressions for G i u i , so that u i is the only variable in the objective. This procedure results in optimal values u 1 * ( m 1 , m 2 ) , u 2 * ( m 1 , m 2 ) for each combination of values for ( m 1 , m 2 ) . The values for x 0 * , x 1 * , x 2 * are determined by these values.
Finally, they iterate over all 1 h 2 of these combinations of values to determine which one maximizes the value of Q . This determines the optimal values of m 1 * , m 2 * for the leader, which in turn determine the optimal strategies for the follower. The overall procedure is summarized in Algorithm 1.
Algorithm 1 Prior approach
for  m 1 = 0 to 1 by increments of h do
       for  m 2 = 0 to 1 by increments of h do
             Calculate values of u 1 , u 2 that minimize | G 1 u 1   | + |   G 2 u 2 | , where 0 u i 1 and the x i satisfy the equations for the equilibrium of the ecological dynamics.
       end for
end for
 Calculate the values for ( m 1 , m 2 ) out of the 1 h 2 possibilities for which the corresponding optimal variables maximize the value of Q .

3. Limitations of Prior Approach

The prior approach has several significant limitations. The first is that it does not check that the calculated values for x i * for given values of u i and m i are biologically sensible. Since the x i correspond to populations, they must be nonnegative.
Another limitation is that the coarse discretization of values for m i means that only a small number of possibilities are considered. This also means that the running time will potentially be large, since we must perform 1 h 2 separate optimizations.
Another significant limitation is that it is assumed that G i is maximized when G i u i = 0 , and the boundary cases when it is maximized at u i = 0 or 1 are ignored.
A final limitation is that the procedure invokes the fminsearch algorithm in Matlab, which is not even guaranteed to find a local minimum, let alone a global minimum.

4. New Approach

We now describe our new approach that addresses the limitations of the prior approach. We will formulate a single quadratic program that corresponds to the full optimization problem and solve it using Gurobi’s nonconvex MIQCP solver which has a guarantee of global optimality (subject to numerical precision).
First we have the main decision variables x i , u i , m i with x i 0 , m i 0 , 0 u i 1 . The objective function Q is a quadratic function of these variables. Next we encode the conditions for equilibrium of the ecological dynamics. We must define several auxiliary variables to do this.
First define η i = g i u i , and τ i = e η i for i = 1 , 2 . For the latter, we use Gurobi’s addGenConstrExp function that uses a piecewise linear approximation for the exponential function. We set these variables to be nonnegative. To provide tighter upper bounds we can set η i g i , τ i e g i , since u i 1 . We next define the auxiliary variable γ i = m i k i + b u i . We can do this by including the quadratic constraint k i γ i + b u i γ i m i = 0 for i = 1 , 2 . Using these variables we can now encode the conditions for equilibrium of ecological dynamics using constraints that are quadratic in the variables.
Next we must encode the conditions that u i is a maximizer of G i . To do this we define several additional auxiliary variables. We define σ i = e g i u i by adding in the constraint σ i τ i = 1 . Next we define β i = 1 ( k i + b i u i ) 2 . We can do this by including the quadratic constraint β i b i 2 u i 2 2 k i b i u i k i 2 = 0 . Finally we define ω i = m i ( k i + b i u i ) 2 by including the quadratic constraint ω i β i m i = 0 . Using these variables, we can now encode the expressions for G i u i that are quadratic in the variables.
Recall that we are trying to select u i [ 0 , 1 ] to maximize G i , for i = 1 , 2 . We can do this by introducing two Lagrange multipliers λ i 1 0 , λ i 2 0 . Then the KKT optimality condition is equivalent to the following three constraints:
G i u i ( u i ) + λ i 1 λ i 2 = 0
λ i 1 ( u i 0 ) = 0
λ i 2 ( u i 1 ) = 0
These constraints are all quadratic in the variables and ensure that we find u i [ 0 , 1 ] that maximizes G i regardless of whether it is at the boundary or at an interior solution with the derivative equal to zero.
Our full formulation is given below. Here the X i j correspond to the same quantities as before and are just defined to simplify presentation, not as new variables.
max m , u , x Q max c ( x 0 2 + x 1 2 + x 2 2 + 2 x 0 x 1 + 2 x 0 x 2 + 2 x 1 x 2 ) K 2 w 1 m 1 2 w 2 m 2 2 r 1 u 1 2 r 2 u 2 2 s . t . x 0 = K ( X 01 + X 02 + X 03 ) r max ( a 0 a 3 ) ( a 0 2 2 a 1 a 2 + a 0 a 3 ) X 01 = a 0 2 a 3 2 d m 1 k 1 m 2 k 2 + r max X 02 = a 1 ( a 3 a 0 ) τ 1 d τ 1 m 2 k 2 + r max τ 1 γ 1 X 03 = a 1 ( a 3 a 0 ) τ 2 d τ 2 m 1 k 1 + r max τ 2 γ 2 x 1 = K ( X 11 + X 12 + X 13 ) r max ( a 0 a 3 ) ( a 0 2 2 a 1 a 2 + a 0 a 3 ) X 11 = a 2 ( a 3 a 0 ) d m 1 k 1 m 2 k 2 + r max X 12 = ( a 0 2 a 1 a 2 ) τ 1 d τ 1 m 2 k 2 + r max τ 1 γ 1 X 13 = ( a 1 a 2 a 0 a 3 ) τ 2 d τ 2 m 1 k 1 + r max τ 2 γ 2 x 2 = K ( X 21 + X 22 + X 23 ) r max ( a 0 a 3 ) ( a 0 2 2 a 1 a 2 + a 0 a 3 ) X 21 = a 2 ( a 3 a 0 ) d m 1 k 1 m 2 k 2 + r max X 22 = ( a 1 a 2 a 0 a 3 ) τ 1 d τ 1 m 2 k 2 + r max τ 1 γ 1 X 23 = ( a 0 2 a 1 a 2 ) τ 2 d τ 2 m 1 k 1 + r max τ 2 γ 2 g i r max σ i 1 α i 0 x 0 + α i 1 x 1 + α i 2 x 2 K + b i ω i + λ i 1 λ i 2 = 0   for   i = 1 , 2 u i λ i 1 = 0   for   i = 1 , 2 u i λ i 2 λ i 2 = 0   for   i = 1 , 2 η i = g i u i   for   i = 1 , 2 τ i = e η i   for   i = 1 , 2 k i γ i + b u i γ i m i = 0   for   i = 1 , 2 σ i τ i = 1   for   i = 1 , 2 β i b i 2 u i 2 2 k i b i u i k i 2 = 0   for   i = 1 , 2 ω i β i m i = 0   for   i = 1 , 2 m i 0   for   i = 1 , 2 x i 0   for   i = 1 , 2 0 u i 1   for   i = 1 , 2 , 3 0 η i g i   for   i = 1 , 2 0 τ i e g i   for   i = 1 , 2 γ i 0   for   i = 1 , 2 σ i 0   for   i = 1 , 2 β i 0   for   i = 1 , 2 ω i 0   for   i = 1 , 2 λ i j 0   for   i = 1 , 2   and   j = 1 , 2
This formulation addresses the limitations of the prior approach. It ensures that all quantities are biologically relevant by imposing nonnegativity constraints on corresponding variables. It allows m i to take on arbitrary nonnegative values, not a small set of discretized values. It involves solving a single optimization problem instead of 1 h 2 separate optimization problems. It uses KKT conditions to ensure that the values of u i that maximize G i are found regardless of whether they are interior or boundary solutions. And the approach guarantees finding a global optimum since that is guaranteed by Gurobi’s nonconvex MIQCP solver.

5. Experiments

We ran experiments with both approaches on a problem instance using the same parameter values as the prior approach [2], which are provided in Table 2. All experiments were done on a single core of a laptop using Windows 11. Experiments with the prior approach used Matlab version 24.1.0.2689473 (R2024a) Update 6 [3], and experiments with the new approach were done using Gurobi version 11.03 [4] with Java version 14.0.2. For the optimizations in the prior approach, Matlab’s function fminsearchbnd was called using parameters TolX = 1 × 10 12 , MaxFunEvals = 1000 . The results are shown in Table 3. We can see that the prior approach found a solution with a negative value for x 1 * , which is not biologically sensible. The prior approach took nearly five minutes while our new approach took less than two seconds.

6. Conclusions

We presented a new approach for computing Stackelberg equilibrium strategies in a Stackelberg evolutionary game dynamic model of cancer evolution previously studied. Our approach is based on solving a new quadratic program formulation. We noted several limitations of the approach used by prior work which are addressed by our approach. When we compared the approaches on a sample instance our approach ran significantly faster and the prior approach output a solution that is not biologically relevant. As more complex game-theoretic and optimization models are being formulated for problems in biology and cancer treatment in particular, it is important to develop efficient algorithms for accurately solving them. While we focused on one instantiation presented in prior work, our approach is applicable more generally to computing optimal strategies in Stackelberg evolutionary games.

Funding

This research received no external funding.

Data Availability Statement

The original contributions presented in the study are included in the article; further inquiries can be directed to the author.

Conflicts of Interest

The author declares no conflicts of interest.

References

  1. Kleshnina, M.; Streipert, S.; Brown, J.; Staňková, K. Game Theory for Managing Evolving Systems: Challenges and Opportunities of Including Vector-Valued Strategies and Life-History Traits. Dyn. Games Appl. 2023, 13, 1130–1155. [Google Scholar] [CrossRef]
  2. Kleshnina, M. Available online: https://github.com/kleshnina/SEGopinion (accessed on 27 August 2024).
  3. The MathWorks, Inc. MATLAB, Version: 24.1.0.2689473 (R2024a) Update 6; The MathWorks, Inc.: Natick, MA, USA, 2024. [Google Scholar]
  4. Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual; Gurobi Optimization LLC: Beaverton, OR, USA, 2024. [Google Scholar]
Table 1. Interpretations of model parameters.
Table 1. Interpretations of model parameters.
ParameterInterpretation
r max Max cell growth rate
g i Cost of resistance strategy (cell type) i
α i j Interaction coefficient between cell types i and j
KCarrying capacity
dNatural death rate
k i Innate resistance that may be present before drug exposure
b i Benefit of the evolved resistance trait in reducing therapy efficacy
Q max Quality of life of a healthy patient
w i Toxicity of drug i
r i Effect of resistance rate of cell type i
cWeight for impact of tumor burden vs. drug toxicities and resistance rates
Table 2. Parameter values used in experiments.
Table 2. Parameter values used in experiments.
ParameterValue
r max 0.45
g 1 0.5
g 2 0.5
a 0 1
a 1 0.15
a 2 0.9
a 3 0.9
K10,000
d0.01
k 1 5
k 2 5
b 1 10
b 2 10
Q max 1
w 1 0.5
w 2 0.2
r 1 0.4
r 2 0.4
c0.5
Table 3. Experimental results for both approaches.
Table 3. Experimental results for both approaches.
Prior ApproachNew Approach
m 1 * 0.40.40837
m 2 * 0.50.46579
u 1 * 0.190150.21361
u 2 * 0.31230.28554
x 0 * 5634.37745749.8474
x 1 * −360.26581.3366
x 2 * 1316.2683950.5000
Q * 0.599360.59780
Running time (seconds)282.691.65
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Ganzfried, S. Computing Stackelberg Equilibrium for Cancer Treatment. Games 2024, 15, 45. https://doi.org/10.3390/g15060045

AMA Style

Ganzfried S. Computing Stackelberg Equilibrium for Cancer Treatment. Games. 2024; 15(6):45. https://doi.org/10.3390/g15060045

Chicago/Turabian Style

Ganzfried, Sam. 2024. "Computing Stackelberg Equilibrium for Cancer Treatment" Games 15, no. 6: 45. https://doi.org/10.3390/g15060045

APA Style

Ganzfried, S. (2024). Computing Stackelberg Equilibrium for Cancer Treatment. Games, 15(6), 45. https://doi.org/10.3390/g15060045

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop