# MetaGPT: Merging Large Language Models
Using Model Exclusive Task Arithmetic

Yuyan Zhou∗
Baichuan Inc.
&Liang Song
Baichuan Inc.
&Bingning Wang
Baichuan Inc. &Weipeng Chen
Baichuan Inc Equal contributionCorresponding author<EMAIL_ADDRESS>
###### Abstract

The advent of large language models (LLMs) like GPT-4 has catalyzed the
exploration of multi-task learning (MTL), in which a single model demonstrates
proficiency across diverse tasks. Task arithmetic has emerged as a cost-
effective approach for MTL. It enables performance enhancement across multiple
tasks by adding their corresponding task vectors to a pre-trained model.
However, the current lack of a method that can simultaneously achieve optimal
performance, computational efficiency, and data privacy limits their
application to LLMs. In this paper, we propose Model Exclusive Task Arithmetic
for merging GPT-scale models (MetaGPT), which formalizes the objective of
model merging into a multi-task learning framework, aiming to minimize the
average loss difference between the merged model and each individual task
model. Since data privacy limits the use of multi-task training data, we
leverage LLMs’ local linearity and task vectors’ orthogonality to separate the
data term and scaling coefficients term and derive a model-exclusive task
arithmetic method.

Our proposed MetaGPT is data-agnostic and bypasses the heavy search process,
making it cost-effective and easy to implement for LLMs. Extensive experiments
demonstrate that MetaGPT leads to improvements in task arithmetic and achieves
state-of-the-art performance on multiple tasks.

##  1 Introduction

Figure 1: Existing methods face the trilemma of performance, data privacy, and
computational costs, which hinders its application to LLMs. Our MetaGPT can
solve these problems under careful approximation and thus can scale to
GPT3-scale LLMs.

In recent years, a well-established paradigm for AI has been to pre-train
models using large-scale datasets and then to fine-tune the models on
different tasks through supervised learning with task-specific datasets, which
can lead to improved performance while requiring less labeled data Devlin et
al. (2018); OpenAI (2023); Dodge et al. (2020); Yang et al. (2023a). However,
for each new application, a separate model has to be fine-tuned and deployed,
which is computationally expensive and resource-intensive Fifty et al. (2021);
Zhang & Yang (2021). Thus, Multi-Task Learning (MTL) methods have been
proposed and developed to enable a single model to solve multiple tasks
concurrently.

Conventional MTL approaches typically involve collecting raw data across
multiple tasks and then jointly training a single model Caruana (1997); Yang
et al. (2023b). However, the fine-tuning process becomes extremely
computationally intensive with the development of large language models (LLMs)
that may comprise billions or even trillions of parameters. Therefore,
researchers have explored merging various task-specific models with the
expectation that the merged model can handle multiple tasks simultaneously.

One of the outstanding merging methods is task arithmetic Ilharco et al.
(2023). For a given task, the element-wise difference between the weights of
the pre-trained model and the fine-tuned model is referred to as the task
vector. Recent studies have shown that linearly adding multiple scaled task
vectors to the pre-trained model can improve performance across those tasks
Ilharco et al. (2023); Yang et al. (2023c). Nevertheless, previous task
arithmetic methods face a trilemma in practice. 1) The best-performing task
arithmetic methods require extra training to obtain optimal hyper-parameters,
but the high computational costs hinder their application to GPT3-scale LLMs.
2) Some training-free methods heuristically set the scaling coefficient to a
constant (e.g., 0.3), which is efficient but leads to sub-optimal performance.
3) Some methods conduct grid search on the training/validation set, which is
sometimes impractical and faces the risk of data privacy concerns. In summary,
as illustrated in Figure 1, there is essentially no task arithmetic method
suitable for billion-scale models that perform satisfactorily in practice.

To address the aforementioned problems, in this paper, we propose MetaGPT: an
optimal and efficient task arithmetic method for MTL without any data (model
exclusive task arithmetic). We begin by providing a detailed theoretical
analysis of the task loss difference and average loss difference introduced by
the task arithmetic algorithm. Since we aim to choose parameters that minimize
the average loss difference, we first separate the data term and scaling
coefficients, which also establishes a performance upper bound for task
arithmetic. After separating the scaling coefficients, the final result is
quadratic for each scaling coefficient, leading to a closed-form solution that
is simple and effective to implement.

The experimental results on the LLaMA-2 Touvron et al. (2023) and Mistral
Jiang et al. (2023) series demonstrate that the MetaGPT approach is superior
to previous merging methods on several tasks. MetaGPT provides an efficient
avenue to optimally implement task arithmetic for large-scale multi-task
learning (MTL) and push the frontiers of language model merging. To sum up,
our contributions include:

  1. 1.

We provide the mathematical formulation of the optimization objective for task
arithmetic and the first theoretical analysis of the performance bound for
task arithmetic.

  2. 2.

To achieve efficient, optimal, and model-exclusive task arithmetic, we
separate the data term and scaling coefficients in the optimization objective,
which leads to a closed-form solution for the scaling coefficients.

  3. 3.

Our MetaGPT is orthogonal to existing task vector-improving methods and can be
integrated to achieve higher performance.

  4. 4.

Extensive experiments demonstrate that our MetaGPT can improve task arithmetic
and achieve state-of-the-art performance.

##  2 Related Work

Model Merging. Currently, model merging has been developed for multiple uses
such as improving performance on a single target task Izmailov et al. (2018);
Wortsman et al. (2022); Zheng et al. (2024), improving out-of-domain
generalization Ramé et al. (2023); Cha et al. (2021); Arpit et al. (2022), and
improving the performance of multi-task learning Ilharco et al. (2023); Yadav
et al. (2024); Yu et al. (2023); Huang et al. (2024); Ye et al. (2023), which
is the core focus of our research. The range of applications has led to a
proliferation of methods to improve beyond simple parameter averaging. Fisher
merging Matena & Raffel (2022) tries to weight the importance of individual
models using Fisher Information Matrix and uses it to merge different models.
RegMean Jin et al. (2022) formulate the merging problem as a regression
problem and leads to an optimal solution for linear models. Task Arithmetic
Ilharco et al. (2023) presents a method for merging models by adding task
vectors to the pre-trained model to improve multi-task performance. Ties
Merging Yadav et al. (2024) and DARE Yu et al. (2023) propose to refine the
task vectors by resolving the interference and removing extremely redundant
components. Ortiz-Jimenez et al. (2024) propose that fine-tuning the models in
their tangent space can amplify weight disentanglement and lead to substantial
performance improvements.

Multi-Task Learning. Multi-task learning is a powerful method for solving
multiple correlated tasks simultaneously Caruana (1997). Current MTL works
mainly focus on learning the shared representations from designing specific
architecture Misra et al. (2016); Sun et al. (2020) or using specific
optimization methods (Sener & Koltun, 2018; Liu et al., 2021). The former
focuses on learning the shared representation using different methods such as
designing specific representation sharing module (Liu et al., 2019; Ding et
al., 2021), learning to branch (Lu et al., 2017; Guo et al., 2020), and based
selection criteria (Ma et al., 2018; Hazimeh et al., 2021). And the latter
focuses on balancing multiple tasks from the perspectives of task training
weights (Sener & Koltun, 2018; Liu et al., 2019), gradient dominance (Chen et
al., 2018; He et al., 2022; Yang et al., 2023b), and solving gradient
conflicts (Yu et al., 2020; Chen et al., 2020; Liu et al., 2021). However, the
conventional MTL approaches for collecting raw data across multiple tasks for
joint training are not suitable for LLMs. The factors contributing to this
issue are twofold: first, computational inefficiency due to the substantial
computational costs associated with updating pre-trained models; second, a
significant number of data proprietors are reluctant to disclose valuable or
privacy-sensitive raw data.

Figure 2: Current task arithmetic based methods face the problems of sub-
optimal performance, huge computational and memory cost, curse of
dimensionality and data privacy, which makes it difficult to scale to LLMs.
Our method solves the aforementioned problems and provides an avenue to scale
task arithmetic to LLMs.

##  3 Preliminaries

###  3.1 Notation

Let $f:\mathcal{X}\times\Theta\to\mathcal{Y}$ be a neural network taking
inputs $\bm{x}\in\mathcal{X}$ and parameterized by a set of weights
$\bm{\theta}\in\Theta$. We assume $\mathcal{X}\subseteq\mathbb{R}^{p}$,
$\Theta\subseteq\mathbb{R}^{m}$ and $\mathcal{Y}\subseteq\mathbb{R}^{q}$. We
consider fine-tuning a pre-trained model $f(\cdot,\bm{\theta}_{0})$ on $T$
different tasks, with each task $t$ consisting of a triplet
$(\mathcal{D}_{t},\mathcal{L}_{t},\bm{\theta}_{t})$ , where
$\mathcal{D}_{t}=(\mathcal{D}^{\text{train}}_{t},\mathcal{D}^{\text{val}}_{t},\mathcal{D}^{\text{test}}_{t})$
is the training, validation and test data of task $t$, $\mathcal{L}_{t}$ is
the loss function of task $t$, and $\bm{\theta}_{t}$ is the model parameters
fine-tuned on task $t$ based on the pre-trained weight $\bm{\theta}_{0}$.

###  3.2 Task Arithmetic

Let the task vector of task $t$ be the difference between the fine-tuned and
the pre-trained weights:

$\bm{\tau}_{t}=\bm{\theta}_{t}-\bm{\theta}_{0}.$ (1)

Task arithmetic aims to solve the multi-task learning problem by directly
adding the scaled task vectors to the pre-trained model weight
$\bm{\theta}_{0}$:

$\bm{\theta}_{\text{final}}=\bm{\theta}_{0}+\sum_{i=1}^{T}\lambda_{i}\bm{\tau}_{i}$
(2)

where $\lambda_{i}$ is the scaling coefficient of task vector $\tau_{i}$. As
illustrated in Eq. 2, the task arithmetic introduces $T$ hyper-parameters
$\\{\lambda_{i}|i=1,\cdots,T\\}$ and the choice of these scaling coefficients
has a significant influence on the performance of the merged model. Thus,
selecting the appropriate scaling coefficients for different task vectors
remains a challenging problem.

###  3.3 Existing Methods

Earlier task arithmetic Ilharco et al. (2023); Yadav et al. (2024) propose to
perform a grid search (G-Task Arithmetic) on the validation set to choose the
optimal scaling coefficients. However, as the number of tasks increases,
exploring all the scaling coefficient combinations faces the curse of
dimensionality. Therefore, to simplify the problem, they use the same value
for multiple scaling coefficients, thereby reducing the computational
complexity. In the absence of the training/validation data, they set
$\lambda=0.3$ as the default setting for dataless arithmetic. Moreover,
Adamerging Yang et al. (2023c) aims to autonomously learn the coefficients
from unlabeled test samples using entropy minimization.

###  3.4 Scalability Challenges for LLMs

The methods mentioned above are not suitable for scaling to LLMs: The grid
search method requires extra validation/training data, which faces the risk of
data privacy concerns and the curse of dimensionality when the number of tasks
increases. For instance, conducting a grid search for three hyper-parameters,
each with a discretization interval of 0.01, would require $10^{6}$ forward
passes across the entire dataset. Setting a fixed value such as $0.3$ for all
the $\lambda_{i}$ is time-efficient and can be applied to LLMs, but it leads
to sub-optimal performance. Using test data input to unsupervised optimize
these hyper-parameters can lead to an optimal solution but requires extra data
and necessitates loading multiple models for training. This process is both
time and memory consuming, making it challenging to apply to LLMs. For
example, merging three LLMs requires loading three LLMs simultaneously to
optimize, which is extremely costly. The statement above suggests that scaling
up existing optimal task arithmetic to LLMs remains a challenging problem.

##  4 Our Proposed MetaGPT

###  4.1 Overview

To solve the problems above, we propose a new algorithm MetaGPT, based on
careful approximations to a closed-form solution, which easily scales to giant
models both in terms of runtime as well as performance while protecting data
privacy. In this section, we state the motivation and optimization problem and
solve it step by step. All proofs of lemmas and theorems are provided in the
appendix.

###  4.2 MetaGPT Optimization Objective

######  Definition 1 (Single Task Loss Difference).

For the fine-tuned model $\bm{\theta}_{i}$ and the task arithmetic merged
model $\bm{\theta}_{\text{final}}$. The Task Loss Difference in task $t$
($\text{TLD}_{t}$) is defined as:

$\displaystyle\textup{TLD}_{t}(\lambda_{1},\cdots,$
$\displaystyle\lambda_{T},\bm{\tau}_{1},\cdots,\bm{\tau}_{T})=\mathcal{L}_{t}(\bm{\theta}_{\textup{final}},\bm{x})-\mathcal{L}_{t}(\bm{\theta}_{t},\bm{x}).$
(3)

It is obvious that smaller $\textup{TLD}_{t}$ suggests that the loss of the
merged model is close or even lower than the fine-tuned model on task $t$,
which indicates a better task arithmetic performance.

However, for task arithmetic, it aims to improve the average performance of
the final model on all the tasks. Thus, we define the average of all the task
loss differences as Average Loss Difference (ALD), which can be formulated as
follows:

######  Definition 2 (Average Task Loss Difference).

For the fine-tuned models $\\{\bm{\theta}_{i}|i=1,\cdots,T\\}$ and task
arithmetic merged model $\bm{\theta}_{\text{final}}$. The average loss
difference for all tasks is defined as:

$\displaystyle\textup{ALD}(\lambda_{1},$
$\displaystyle\cdots,\lambda_{T},\bm{\tau}_{1},\cdots,\bm{\tau}_{T})=\frac{1}{T}\sum_{t=1}^{T}\left(\mathcal{L}_{t}(\bm{\theta}_{\textup{final}},\bm{x})-\mathcal{L}_{t}(\bm{\theta}_{t},\bm{x})\right).$
(4)

Thus, the optimization objective of MetaGPT is to find the optimal scaling
coefficients that can minimize the ALD, which can be formulated as:

######  Definition 3 (Optimization objective of MetaGPT).

Our MetaGPT aims at finding the scaling coefficients
$\\{\lambda_{i}|i=1,\cdots,T\\}$, which minimizes the average loss difference
ALD:

$\displaystyle\mathop{\arg\min}_{\lambda_{1},\cdots,\lambda_{T}}\frac{1}{T}\sum_{t=1}^{T}\left(\mathcal{L}_{t}(\bm{\theta}_{\textup{final}},\bm{x})-\mathcal{L}_{t}(\bm{\theta}_{t},\bm{x})\right).$
(5)

###  4.3 Separating Data and Coefficients

Before analyzing ALD, we start with reformulating $\textup{TLD}_{t}$ by its
Taylor expansion.

######  Lemma 4.

Using Taylor expansion for $\mathcal{L}(\bm{\theta}_{\textup{final}},\bm{x})$
at $\bm{\theta}_{t}$, the $\textup{TLD}_{t}$ in Eq. 3 can be reformulated as a
quadratic form with respect to the linear combination of $\bm{\lambda}$ and
$\bm{\theta}$:

$\displaystyle\textup{TLD}_{t}$
$\displaystyle=\frac{1}{2}\bm{h}_{t}^{\top}\left(\int_{0}^{1}\nabla^{2}\mathcal{L}_{t}(\gamma_{t}(\beta))d\beta\right)\bm{h}_{t},$
(6)

where
$\gamma_{t}(\beta)=\bm{\theta}_{t}+\beta(\bm{\theta}_{\textup{final}}-\bm{\theta}_{t})$
and $\bm{h}_{t}$ is the linear combination of $\bm{\lambda}$ and
$\bm{\theta}$:

$\displaystyle\bm{h}_{t}=\sum_{k\neq
t}\lambda_{k}(\bm{\theta}_{k}-\bm{\theta}_{0})-(1-\lambda_{t})(\bm{\theta}_{t}-\bm{\theta}_{0}).$
(7)

Single $\textup{TLD}_{t}$ is associated with the data, models, and scaling
coefficients. As we can see in Eq. 6, we have transformed the data term
$\bm{x}_{t}$ to the Hessian, the coefficients
$\bm{\lambda}=[\lambda_{1},\cdots,\lambda_{T}]$ and models term
$[\bm{\theta}_{1},\cdots,\bm{\theta}_{T}]$ to $\bm{h}$. As our method tends to
achieve model-exclusive task arithmetic, the final result should not correlate
with the data term. Thus, we first provide a property, which will be used
latter in our theorem proofs to separate the data term and scaling
coefficients and models term. In general, if a pre-trained network
$f(\cdot;\bm{\theta}_{0})$ demonstrates kernel behavior during fine-tuning,
i.e., fine-tuning occurs in the linear regime, the following property must be
satisfied Jacot et al. (2018):

######  Property 5 (NTK linearization).

Around the initialization weights $\bm{\theta}_{0}$, a neural network can be
approximated with a linear approximation:

$\displaystyle
f(\bm{x};\bm{\theta}_{0}+\alpha(\bm{\theta}_{t}-\bm{\theta}_{0}))\approx
f(\bm{x};\bm{\theta}_{0})+\alpha\cdot C.$ (8)

where $C=(\bm{\theta}_{t}-\bm{\theta}_{0})^{\top}\nabla
f(\bm{x},\bm{\theta}_{0})$ is a data and model dependent constant.

It is worth noting that, as the network width approaches infinity, Eq. 8
becomes exact and remains valid throughout training Jacot et al. (2018); Arora
et al. (2019); Lee et al. (2019), which is specifically suitable for the LLMs
arithmetic scenario.

The second property is observed by (Ilharco et al., 2023), which states that
the different task vectors are orthogonal:

######  Property 6 (Orthogonality of Task Vectors).

For task vector $\bm{\tau}_{i}=\bm{\theta}_{i}-\bm{\theta}_{0}$ and
$\bm{\tau}_{j}=\bm{\theta}_{j}-\bm{\theta}_{0}$ ($i\neq j$), we have the
following equation:

$\bm{\tau}_{i}^{\top}\bm{\tau}_{j}=(\bm{\theta}_{i}-\bm{\theta}_{0})^{\top}(\bm{\theta}_{j}-\bm{\theta}_{0})=0.$
(9)

Now, as we previously introduce our first Lemma to transform the
$\textup{TLD}_{t}$ in Eq. 3 into a quadratic form with respect to the linear
combination of $\bm{\lambda}$ and $\bm{\theta}$. Next, using Property 5,6 and
Lemma 7, we can upper bound the $\textup{TLD}_{t}$ and separate the data term
and scaling coefficients and models term.

######  Theorem 7.

The $\textup{TLD}_{t}$ can be upper bounded by:

$\displaystyle\textup{TLD}_{t}(\lambda_{1},\cdots,\lambda_{T},\bm{\tau}_{1},\cdots,\bm{\tau}_{T})\leq\frac{\delta_{t}^{2}}{2}\left\|\bm{\theta}_{t}-\bm{\theta}_{0}\right\|_{2}^{2}\bigg{\\{}\sum_{k\neq
t}^{T}\vmathbb{1}_{t}(\lambda_{k}^{2})\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\bigg{\\}}$

where $\delta_{t}$ is a data-dependent constant and we use
$\vmathbb{1}_{t}(\lambda_{k}^{2})$ to denote
$(\lambda^{2}_{k})\vmathbb{1}(k\neq t)+(1-\lambda^{2}_{k})\vmathbb{1}(k=t)$.
Now, after separating the data-related term to $\delta_{t}$, the scaling
coefficients and models term to $\vmathbb{1}_{t}(\lambda_{k}^{2})$. By summing
all the $\textup{TLD}_{t}$s, we can separate the two terms for ALD:

######  Theorem 8.

By summing all the $\textup{TLD}_{t}$, we can separate the correlation between
data term and scaling coefficients term in ALD:

ALD
$\displaystyle(\lambda_{1},\cdots,\lambda_{T},\bm{\tau}_{1},\cdots,\bm{\tau}_{T})\leq\sum_{t=1}^{T}\delta_{t}^{2}\left\|\bm{\theta}_{t}-\bm{\theta}_{0}\right\|_{2}^{2}\left\\{\sum_{k\neq
t}^{T}\vmathbb{1}(\lambda_{k}^{2})\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\right\\},$
(10)

###  4.4 The Optimal Solution

After separating the data term and the scaling coefficients term, we can now
reformulate our optimization objective Eq. 10 and derive the closed-form
optimal solution of the scaling coefficients.

######  Theorem 9 ($\lambda$ decomposition of ALD).

For each $\lambda_{t}$, we use it to decompose Eq. 10 as:

$\displaystyle\textup{ALD}\leq\sum_{t=1}^{T}\textup{ALD}_{\lambda_{t}},$ (11)

where $\textup{ALD}_{\lambda_{t}}$ is:

$\displaystyle\textup{ALD}_{\lambda_{t}}=$
$\displaystyle\frac{\delta_{0}^{2}}{2}\|\bm{\theta}_{t}-\bm{\theta}_{0}\|^{2}\left[\sum_{k=1}^{T}\vmathbb{1}_{t}(\lambda)\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\right],$
(12)

where $\delta_{0}=\max_{t}\delta_{t}$. The equation above easily leads to a
model-exclusive closed-form solution:

######  Theorem 10 (Optimal Scaling Coefficients).

We can solve $\lambda_{t}$ form Eq 12 by:

$\displaystyle\lambda_{t}=\mathop{\arg\min}_{\lambda_{t}}\|\bm{\theta}_{t}-\bm{\theta}_{0}\|^{2}\left[\sum_{k=1}^{T}\vmathbb{1}_{t}(\lambda)\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\right].$
(13)

The above equation is quadratic on $\lambda_{t}$ and the optimal solution for
$\lambda_{t}$ is:

$\displaystyle\lambda_{t}=\frac{\|\bm{\theta}_{t}-\bm{\theta}_{0}\|^{2}}{\sum_{k=1}^{n}\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}}.$
(14)

Figure 3: Verification of NTK linearization. We randomly sampled the outputs
of Llama-2-7b-chat-hf with different $\alpha$. We can see that the sampled
outputs are linearly with $\alpha$ as expected.

Figure 4: Verification of orthogonality. We calculate the cosine similarity
between six different task vectors and find that their cosine similarity is
nearly 0.

##  5 Property Verification

In Section 4, we introduced two properties essential to our proof. In this
section, we conduct experiments to verify these properties.

###  5.1 NTK Linearization

Jacot et al. (2018) have proved that when the width of the neural network
approaches infinity, it demonstrates kernel behavior and the optimization
proceeds in the linear regime. We test Llama-2-7b-chat-hf Touvron et al.
(2023) on AGIEval Zhong et al. (2023) dataset to verify its linearity. We have
randomly sampled three outputs of the Llama-2-7b-chat-hf when $\alpha$ in Eq.
8 gets value of $[0,0.1,\cdots,1]$. For better visualization, we also subtract
all the outputs using $\max\\{y_{i}\\}$, ensuring they have the same endpoint.
From the results in Figure 4, we can see that all the outputs are almost
linear with $\alpha$, which indicates that LLMs do exhibit a kernel behavior
during finetuning.

###  5.2 Task Vector Orthogonality

Ilharco et al. (2023); Yang et al. (2023c) have performed experiments to
verify this property for vision models. For LLMs, we also observe similar
results: these task vectors are almost orthogonal to each other. The result
has been shown in Figure 4. We can see that different task vectors are almost
orthogonal, and their cosine similarity is nearly 0 as Eq.9 expected, which
verifies the property we have used for our proof.

Table 1:  Performance comparison of merging different LLaMA-2-7B fine-tuned
models on different datasets.

Model | WinoGrande | AGIEval | GSM8k | MATH | MBPP | HumanEval | Abs. Avg | Nor. Avg
---|---|---|---|---|---|---|---|---
LM | 62.67 | 34.01 | 28.66 | 4.00 | 22.00 | 7.31 | 26.44 | 0.91
Math | 61.64 | 29.40 | 47.16 | 2.40 | 17.40 | 11.58 | 28.26 | 0.84
Code | 61.88 | 27.41 | 17.21 | 2.20 | 24.80 | 21.92 | 25.90 | 0.84
Weight Average | 63.93 | 31.36 | 37.68 | 7.00 | 23.40 | 20.12 | 30.58 | 1.25
Task Arithmetic | 63.54 | 31.70 | 37.53 | 5.20 | 23.20 | 19.51 | 30.11 | 1.12
Ties Merging | 62.67 | 32.10 | 37.93 | 7.40 | 22.80 | 18.29 | 30.20 | 1.26
DARE | 63.27 | 32.25 | 37.86 | 7.00 | 24.40 | 19.51 | 30.72 | 1.26
MetaGPT(ours) | 64.25 | 32.71 | 45.41 | 7.80 | 21.20 | 17.68 | 31.51 | 1.31

Table 2: Performance comparison of merging different Mistral-7B fine-tuned
models on different datasets.

Model | WinoGrande | AGIEval | GSM8k | MATH | MBPP | HumanEval | Abs. Avg | Nor. Avg
---|---|---|---|---|---|---|---|---
LM | 69.30 | 37.55 | 47.54 | 7.80 | 34.40 | 34.75 | 38.56 | 0.776
Math | 63.46 | 38.06 | 68.46 | 28.00 | 24.00 | 25.00 | 41.16 | 0.854
Code | 67.32 | 40.69 | 60.73 | 15.60 | 43.40 | 39.02 | 44.46 | 0.917
Weight Average | 67.88 | 41.12 | 62.77 | 17.40 | 40.20 | 38.41 | 44.63 | 0.921
Task Arithmetic | 67.88 | 41.41 | 63.38 | 18.80 | 40.20 | 38.40 | 45.01 | 0.932
Ties Merging | 67.72 | 41.06 | 60.35 | 17.80 | 40.20 | 40.24 | 44.56 | 0.924
DARE | 67.40 | 40.58 | 59.67 | 19.00 | 36.00 | 40.85 | 43.92 | 0.913
MetaGPT(ours) | 68.35 | 41.86 | 66.03 | 20.80 | 39.00 | 35.37 | 45.24 | 0.936

##  6 Experiments

In this section, we conduct experiments to demonstrate the effectiveness of
our MetaGPT. In the first section, we demonstrate that our MetaGPT
consistently achieves optimal average performance across diverse datasets and
is robust for model series with varying parameter sizes and architectures.
DARE and Ties-Merging are task vector-improving methods that resolve conflicts
and redundant parameters between task vectors. We conduct experiments to
demonstrate that our method is orthogonal to theirs and can be integrated to
improve the average performance further. Finally, we show that the model
merged by our MetaGPT has better out-of-distribution generalization ability.

Table 3: Comparison of performance of merging fine-tuned LLaMA-2-13B on
different datasets.

Model | WinoGrande | AGIEval | GSM8K | MATH | MBPP | HumanEval | Abs. Avg | Nor. Avg
---|---|---|---|---|---|---|---|---
LM | 64.80 | 35.04 | 42.84 | 4.80 | 27.00 | 15.24 | 31.62 | 1.02
Math | 60.38 | 36.74 | 55.27 | 3.40 | 22.60 | 12.80 | 31.87 | 0.93
Code | 63.93 | 32.04 | 36.47 | 5.00 | 26.60 | 16.46 | 30.08 | 1.01
Weight Average | 64.88 | 37.23 | 53.15 | 7.60 | 29.80 | 21.95 | 35.77 | 1.29
Task Arithmetic | 65.11 | 35.48 | 50.34 | 7.20 | 29.80 | 21.95 | 34.98 | 1.25
Ties Merging | 65.23 | 36.02 | 51.23 | 7.40 | 30.20 | 23.17 | 35.54 | 1.28
DAREs | 65.70 | 36.87 | 51.85 | 7.60 | 30.00 | 22.56 | 35.76 | 1.29
MetaGPT(ours) | 65.04 | 37.33 | 52.92 | 7.80 | 30.40 | 21.95 | 35.91 | 1.30

Table 4: MetaGPT can be integrated with DARE and Ties-Merging, thereby leading
to further improvment.

Method | WinoGrande | AGIEval | GSM8k | MATH | MBPP | HumanEval | Abs. Avg | Nor. Avg
---|---|---|---|---|---|---|---|---
Ties-Merging | 62.67 | 32.10 | 37.93 | 7.40 | 22.80 | 18.29 | 30.20 | 1.26
Ties + MetaGPT | 62.35 | 32.91 | 46.10 | 8.00 | 22.40 | 17.68 | 31.57 | 1.33
Dare | 63.27 | 32.25 | 37.86 | 7.00 | 24.40 | 19.51 | 30.72 | 1.26
Dare + MetaGPT | 62.99 | 33.01 | 45.72 | 7.60 | 21.80 | 18.29 | 31.57 | 1.30

###  6.1 Merging Models Using MetaGPT

#### Dataset and Models.

To test the effectiveness of our method, we use Llama-2-7b-chat-hf Touvron et
al. (2023), MAmmoTH-7B Yue et al. (2023) and llama-2-coder-7b Manuel Romero
(2023) as models fine-tuned on general knowledge, math, and code datasets
using the pre-trained model Llama-2-7B-hf Touvron et al. (2023). Moreover, we
use a different model architecture: Mistral-7B-Instruct-v0.2 AI ,
MAmmoTH2-7B-Plus Yue et al. (2024) and Mistral-7B-codealpaca-lora Nondzu as
models fine-tuned on general knowledge, math, and code datasets using pre-
trained model Mistral 7B Jiang et al. (2023). We also provide experiments
using models with larger sizes: Llama-2-13b-chat-hf Touvron et al. (2023),
MAmmoTH-13B Yue et al. (2023), and llama-2-13b-code-chat TAŞAR (2023) as
models fine-tuned on general knowledge, math, and code datasets using the pre-
trained model Llama-2-13B-hf Touvron et al. (2023). We use WinoGrande
Sakaguchi et al. (2021) and AGIEval Zhong et al. (2023) for evaluating general
knowledge performance, GSM8K Cobbe et al. (2021) and MATH Saxton & Hill (2019)
for testing mathematical reasoning ability, HumanEval Chen et al. (2021) and
MBPP Austin et al. (2021) for estimating code-generation capacity.

#### Evaluation Metrics.

We use common evaluation settings for a single task: 5-shot accuracy for
AGIEval, 4-shot accuracy for GSM8K and MATH, 3-shot accuracy for MBPP, and
zero-shot accuracy for HumanEval and WinoGrande. We employ two key metrics in
evaluating different merging methods: absolute average performance and
normalized average accuracy.

#### Quantitative Evaluation for LLaMA-2-7B.

We use the metrics and datasets we introduced above to evaluate the
performance of different methods. We use Weight Average Wortsman et al.
(2022), Task Arithmetic Ilharco et al. (2023), Ties-Merging Yadav et al.
(2024) and DARE Yu et al. (2023), which are also model exclusive and
computationally efficient methods, to compare with our method by merging
LLaMA-2-7B. The scores in Table 1 show that for WinoGrande, AGIEval, GSM8k,
and MATH dataset, our method scores 64.25, 32.71, 45.41, and 7.80, which
outperforms other methods. For the HumanEval dataset, DARE performs best, and
for the MBPP dataset, the Weight Average method achieves the highest score.
Since our method aims to achieve the average best performance, we use absolute
average performance score and normalized average performance score to compare
the five methods. We can see that our MetaGPT achieves the rank-1 score 31.51,
1.31 in both absolute average performance and normalized average performance.

#### Using Different Model Architecture.

We also use a different model architecture, Mistral-7B, for evaluation, and
the result has been shown in Table 2. The scores in Table 2 show similar
results to LLaMA-2-7B: For WinoGrande, AGIEval, GSM8k, and MATH dataset, our
MetaGPT scores 41.86, 68.35, 66.03, 20.8, which outperforms existing methods,
for HumanEval dataset Weight Average, Task Arithmetic, and Ties Merging
performs best and for MBPP dataset, DARE method achieves the highest score.

#### Using Larger Model Size.

We also test our method using a larger model LLaMA-2-13B Touvron et al.
(2023). The scores in Table 3 demonstrate that for AGIEval, Math, and MBPP
datasets, our method outperforms other methods. For WinoGrand, GSM8K, and
HumanEval dataset, DARE, Weight Average and Ties-Merging achieves the highest
score. Similarly, under the average measure absolute average performance and
normalized average performance, our method also outperforms the other five
methods.

#### Integrate with Ties/DARE

As there are conflicts and redundant parameters between task vectors, DARE Yu
et al. (2023) and Ties-Merging Yadav et al. (2024) are two methods trying to
solve the interfaces, reducing the redundancy and thereby improving the
performance of task arithmetic. Since our method is also based on the
framework of task arithmetic, Ties-merging and DARE are expected to improve
the performance of our MetaGPT further. As we can see in Table 4, under the
baseline of Ties-Merging and DARE methods, our method is orthogonal to Ties-
Merging and DARE and can integrate them into our MetaGPT, thus leading to
further improvement. For example, the average absolute performance of DARE has
been improved by our MetaGPT from 30.72 to 31.57. And the normalized absolute
performance of DARE has been improved by our MetaGPT from 1.26 to 1.3. Ties-
merging also leads to a similar conclusion: the average absolute performance
of DARE has been improved by our MetaGPT from 30.20 to 31.57. And the
normalized absolute performance of DARE has been improved by our MetaGPT from
1.26 to 1.33.

###  6.2 Out of Distribution Generalization

Following Yang et al. (2023c); Jin et al. (2022), we also compare the out-of-
distribution generalization ability of different merging methods. We evaluate
different methods using JEC-QA Zhong et al. (2020), FinanceIQ DI (2023), and
MedQA Jin et al. (2021) dataset. All three datasets use 5-shot accuracy as the
evaluation metric. Table 5 summarizes out-of-distribution generalization
performance when merging all domain specific models using different methods.
As we can see, MetaGPT outperforms current methods on these unseen datasets,
which demonstrates that MetaGPT is more robust to the test data distribution
shifts.

Table 5: Out of Distribution Generalization

Model |  JEC-QA |  FinancelQ |  MedQA |  Avg
---|---|---|---|---
LM |  31.32 |  32.83 |  30.20 |  31.45
Math |  25.56 |  30.25 |  24.73 |  26.85
Code |  29.23 |  30.87 |  26.25 |  28.78
Weight Average |  30.73 |  34.17 |  29.90 |  31.60
Task Arithmetic |  30.85 |  33.89 |  30.13 |  31.62
Ties Merging |  30.80 |  33.53 |  30.02 |  31.45
DARE |  30.79 |  33.93 |  30.17 |  31.63
MetaGPT(ours) |  30.97 |  34.31 |  30.07 |  31.78

##  7 Conclusion

In this paper, we have provided a novel model merging method named MetaGPT, an
efficient and optimal model-exclusive task arithmetic specifically designed
for LLMs. We provide the mathematical formulation of task arithmetic’s
optimization objective and the theoretical analysis of the task arithmetic
performance bound. By separating the data and scaling coefficient term under
careful approximation, the closed-form solution provides an avenue for
optimally achieving task arithmetic without using any data. Extensive
experiment results show that our MetaGPT outperforms the existing state-of-
the-art model-exclusive merging method and can be integrated with task vector-
improving methods such as Ties-Merging and DARE.

##  8 Limitations

(1) Our works share the same general limitation of existing task arithmetic
based methods: Our merging method relies on common initialization and model
architecture, which ensures that the task vectors are orthogonal. (2)
Moreover, since our method is specifically designed for LLMs and relies on the
NTK linearization, for small size models, our method may not perform well.

## References

  * (1) Mistral AI.  Mistral-7b-instruct-v0.2.  URL https://huggingface.co/mistralai/Mistral-7B-Instruct-v0.2.
  * Arora et al. (2019) Sanjeev Arora, Simon S. Du, Wei Hu, Zhiyuan Li, Ruslan Salakhutdinov, and Ruosong Wang.  On exact computation with an infinitely wide neural net.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019.
  * Arpit et al. (2022) Devansh Arpit, Huan Wang, Yingbo Zhou, and Caiming Xiong.  Ensemble of averages: Improving model selection and boosting performance in domain generalization.  _Advances in Neural Information Processing Systems_ , 35:8265–8277, 2022.
  * Austin et al. (2021) Jacob Austin, Augustus Odena, Maxwell Nye, Maarten Bosma, Henryk Michalewski, David Dohan, Ellen Jiang, Carrie Cai, Michael Terry, Quoc Le, et al.  Program synthesis with large language models.  _arXiv preprint arXiv:2108.07732_ , 2021.
  * Caruana (1997) Rich Caruana.  Multitask learning.  _Machine learning_ , 28:41–75, 1997.
  * Cha et al. (2021) Junbum Cha, Sanghyuk Chun, Kyungjae Lee, Han-Cheol Cho, Seunghyun Park, Yunsung Lee, and Sungrae Park.  Swad: Domain generalization by seeking flat minima.  _Advances in Neural Information Processing Systems_ , 34:22405–22418, 2021.
  * Chen et al. (2021) Mark Chen, Jerry Tworek, Heewoo Jun, Qiming Yuan, Henrique Ponde de Oliveira Pinto, Jared Kaplan, Harri Edwards, Yuri Burda, Nicholas Joseph, Greg Brockman, Alex Ray, Raul Puri, Gretchen Krueger, Michael Petrov, Heidy Khlaaf, Girish Sastry, Pamela Mishkin, Brooke Chan, Scott Gray, Nick Ryder, Mikhail Pavlov, Alethea Power, Lukasz Kaiser, Mohammad Bavarian, Clemens Winter, Philippe Tillet, Felipe Petroski Such, Dave Cummings, Matthias Plappert, Fotios Chantzis, Elizabeth Barnes, Ariel Herbert-Voss, William Hebgen Guss, Alex Nichol, Alex Paino, Nikolas Tezak, Jie Tang, Igor Babuschkin, Suchir Balaji, Shantanu Jain, William Saunders, Christopher Hesse, Andrew N. Carr, Jan Leike, Josh Achiam, Vedant Misra, Evan Morikawa, Alec Radford, Matthew Knight, Miles Brundage, Mira Murati, Katie Mayer, Peter Welinder, Bob McGrew, Dario Amodei, Sam McCandlish, Ilya Sutskever, and Wojciech Zaremba.  Evaluating large language models trained on code, 2021.
  * Chen et al. (2018) Zhao Chen, Vijay Badrinarayanan, Chen-Yu Lee, and Andrew Rabinovich.  Gradnorm: Gradient normalization for adaptive loss balancing in deep multitask networks.  In _ICML_ , pp. 794–803. PMLR, 2018.
  * Chen et al. (2020) Zhao Chen, Jiquan Ngiam, Yanping Huang, Thang Luong, Henrik Kretzschmar, Yuning Chai, and Dragomir Anguelov.  Just pick a sign: Optimizing deep multitask models with gradient sign dropout.  In _NeurIPS_ , 2020.
  * Cobbe et al. (2021) Karl Cobbe, Vineet Kosaraju, Mohammad Bavarian, Mark Chen, Heewoo Jun, Lukasz Kaiser, Matthias Plappert, Jerry Tworek, Jacob Hilton, Reiichiro Nakano, Christopher Hesse, and John Schulman.  Training verifiers to solve math word problems.  _arXiv preprint arXiv:2110.14168_ , 2021.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_ , 2018.
  * DI (2023) Duxiaoman DI.  Financeiq, 2023.  URL https://huggingface.co/datasets/Duxiaoman-DI/FinanceIQ.
  * Ding et al. (2021) Ke Ding, Xin Dong, Yong He, Lei Cheng, Chilin Fu, Zhaoxin Huan, Hai Li, Tan Yan, Liang Zhang, Xiaolu Zhang, et al.  Mssm: a multiple-level sparse sharing model for efficient multi-task learning.  In _SIGIR_ , pp. 2237–2241, 2021.
  * Dodge et al. (2020) Jesse Dodge, Gabriel Ilharco, Roy Schwartz, Ali Farhadi, Hannaneh Hajishirzi, and Noah Smith.  Fine-tuning pretrained language models: Weight initializations, data orders, and early stopping.  _arXiv preprint arXiv:2002.06305_ , 2020.
  * Fifty et al. (2021) Chris Fifty, Ehsan Amid, Zhe Zhao, Tianhe Yu, Rohan Anil, and Chelsea Finn.  Efficiently identifying task groupings for multi-task learning.  _Advances in Neural Information Processing Systems_ , 34:27503–27516, 2021.
  * Guo et al. (2020) Pengsheng Guo, Chen-Yu Lee, and Daniel Ulbricht.  Learning to branch for multi-task learning.  In _ICML_ , pp. 3854–3863. PMLR, 2020.
  * Hazimeh et al. (2021) Hussein Hazimeh, Zhe Zhao, Aakanksha Chowdhery, Maheswaran Sathiamoorthy, Yihua Chen, Rahul Mazumder, Lichan Hong, and Ed Chi.  Dselect-k: Differentiable selection in the mixture of experts with applications to multi-task learning.  _NeurIPS_ , 34:29335–29347, 2021.
  * He et al. (2022) Yun He, Xue Feng, Cheng Cheng, Geng Ji, Yunsong Guo, and James Caverlee.  Metabalance: Improving multi-task recommendations via adapting gradient magnitudes of auxiliary tasks.  _WWW_ , pp. 2205–2215, 2022.
  * Huang et al. (2024) Chenyu Huang, Peng Ye, Tao Chen, Tong He, Xiangyu Yue, and Wanli Ouyang.  Emr-merging: Tuning-free high-performance model merging.  _arXiv preprint arXiv:2405.17461_ , 2024.
  * Ilharco et al. (2023) Gabriel Ilharco, Marco Tulio Ribeiro, Mitchell Wortsman, Suchin Gururangan, Ludwig Schmidt, Hannaneh Hajishirzi, and Ali Farhadi.  Editing models with task arithmetic.  _The Twelfth International Conference on Learning Representations_ , 2023.
  * Izmailov et al. (2018) Pavel Izmailov, Dmitrii Podoprikhin, Timur Garipov, Dmitry Vetrov, and Andrew Gordon Wilson.  Averaging weights leads to wider optima and better generalization.  _arXiv preprint arXiv:1803.05407_ , 2018.
  * Jacot et al. (2018) Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  _Advances in neural information processing systems_ , 31, 2018.
  * Jiang et al. (2023) Albert Q Jiang, Alexandre Sablayrolles, Arthur Mensch, Chris Bamford, Devendra Singh Chaplot, Diego de las Casas, Florian Bressand, Gianna Lengyel, Guillaume Lample, Lucile Saulnier, et al.  Mistral 7b.  _arXiv preprint arXiv:2310.06825_ , 2023.
  * Jin et al. (2021) Di Jin, Eileen Pan, Nassim Oufattole, Wei-Hung Weng, Hanyi Fang, and Peter Szolovits.  What disease does this patient have? a large-scale open domain question answering dataset from medical exams.  _Applied Sciences_ , 11(14):6421, 2021.
  * Jin et al. (2022) Xisen Jin, Xiang Ren, Daniel Preotiuc-Pietro, and Pengxiang Cheng.  Dataless knowledge fusion by merging weights of language models.  _arXiv preprint arXiv:2212.09849_ , 2022.
  * Lee et al. (2019) Jaehoon Lee, Lechao Xiao, Samuel Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl-Dickstein, and Jeffrey Pennington.  Wide neural networks of any depth evolve as linear models under gradient descent.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019.
  * Liu et al. (2021) Bo Liu, Xingchao Liu, Xiaojie Jin, Peter Stone, and Qiang Liu.  Conflict-averse gradient descent for multi-task learning.  _NeurIPS_ , 34:18878–18890, 2021.
  * Liu et al. (2019) Shikun Liu, Edward Johns, and Andrew J. Davison.  End-to-end multi-task learning with attention.  In _CVPR_ , pp. 1871–1880. Computer Vision Foundation / IEEE, 2019.
  * Lu et al. (2017) Yongxi Lu, Abhishek Kumar, Shuangfei Zhai, Yu Cheng, Tara Javidi, and Rogerio Feris.  Fully-adaptive feature sharing in multi-task networks with applications in person attribute classification.  In _CVPR_ , pp. 5334–5343, 2017.
  * Ma et al. (2018) Jiaqi Ma, Zhe Zhao, Xinyang Yi, Jilin Chen, Lichan Hong, and Ed H. Chi.  Modeling task relationships in multi-task learning with multi-gate mixture-of-experts.  In _SIGKDD_ , pp. 1930–1939. ACM, 2018.
  * Manuel Romero (2023) Manuel Romero.  llama-2-coder-7b (revision d30d193), 2023.  URL https://huggingface.co/mrm8488/llama-2-coder-7b.
  * Matena & Raffel (2022) Michael S Matena and Colin A Raffel.  Merging models with fisher-weighted averaging.  _Advances in Neural Information Processing Systems_ , 35:17703–17716, 2022.
  * Misra et al. (2016) Ishan Misra, Abhinav Shrivastava, Abhinav Gupta, and Martial Hebert.  Cross-stitch networks for multi-task learning.  In _CVPR_ , pp. 3994–4003. IEEE Computer Society, 2016.
  * (34) Nondzu.  Mistral-7b-codealpaca-lora.  URL https://huggingface.co/Nondzu/Mistral-7B-codealpaca-lora.
  * OpenAI (2023) OpenAI.  GPT-4 technical report, 2023.
  * Ortiz-Jimenez et al. (2024) Guillermo Ortiz-Jimenez, Alessandro Favero, and Pascal Frossard.  Task arithmetic in the tangent space: Improved editing of pre-trained models.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Ramé et al. (2023) Alexandre Ramé, Kartik Ahuja, Jianyu Zhang, Matthieu Cord, Léon Bottou, and David Lopez-Paz.  Model ratatouille: Recycling diverse models for out-of-distribution generalization.  In _International Conference on Machine Learning_ , pp. 28656–28679. PMLR, 2023.
  * Sakaguchi et al. (2021) Keisuke Sakaguchi, Ronan Le Bras, Chandra Bhagavatula, and Yejin Choi.  Winogrande: An adversarial winograd schema challenge at scale.  _Communications of the ACM_ , 64(9):99–106, 2021\.
  * Saxton & Hill (2019) Grefenstette Saxton and Kohli Hill.  Analysing mathematical reasoning abilities of neural models.  _arXiv:1904.01557_ , 2019.
  * Sener & Koltun (2018) Ozan Sener and Vladlen Koltun.  Multi-task learning as multi-objective optimization.  In _NeurIPS_ , pp. 525–536, 2018.
  * Sun et al. (2020) Ximeng Sun, Rameswar Panda, Rogerio Feris, and Kate Saenko.  Adashare: Learning what to share for efficient deep multi-task learning.  _NeurIPS_ , 33:8728–8740, 2020.
  * TAŞAR (2023) Davut Emre TAŞAR.  llama-2-13b-code-chat, 2023.  URL https://huggingface.co/emre/llama-2-13b-code-chat.
  * Touvron et al. (2023) Hugo Touvron, Louis Martin, Kevin Stone, Peter Albert, Amjad Almahairi, Yasmine Babaei, Nikolay Bashlykov, Soumya Batra, Prajjwal Bhargava, Shruti Bhosale, et al.  Llama 2: Open foundation and fine-tuned chat models.  _arXiv preprint arXiv:2307.09288_ , 2023.
  * Wortsman et al. (2022) Mitchell Wortsman, Gabriel Ilharco, Samir Ya Gadre, Rebecca Roelofs, Raphael Gontijo-Lopes, Ari S Morcos, Hongseok Namkoong, Ali Farhadi, Yair Carmon, Simon Kornblith, et al.  Model soups: averaging weights of multiple fine-tuned models improves accuracy without increasing inference time.  In _International conference on machine learning_ , pp. 23965–23998. PMLR, 2022.
  * Yadav et al. (2024) Prateek Yadav, Derek Tam, Leshem Choshen, Colin A Raffel, and Mohit Bansal.  Ties-merging: Resolving interference when merging models.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Yang et al. (2023a) Aiyuan Yang, Bin Xiao, Bingning Wang, Borong Zhang, Ce Bian, Chao Yin, Chenxu Lv, Da Pan, Dian Wang, Dong Yan, et al.  Baichuan 2: Open large-scale language models.  _arXiv preprint arXiv:2309.10305_ , 2023a.
  * Yang et al. (2023b) Enneng Yang, Junwei Pan, Ximei Wang, Haibin Yu, Li Shen, Xihua Chen, Lei Xiao, Jie Jiang, and Guibing Guo.  Adatask: A task-aware adaptive learning rate approach to multi-task learning.  In _AAAI_ , volume 37, pp. 10745–10753, 2023b.
  * Yang et al. (2023c) Enneng Yang, Zhenyi Wang, Li Shen, Shiwei Liu, Guibing Guo, Xingwei Wang, and Dacheng Tao.  Adamerging: Adaptive model merging for multi-task learning.  In _The Twelfth International Conference on Learning Representations_ , 2023c.
  * Ye et al. (2023) Peng Ye, Chenyu Huang, Mingzhu Shen, Tao Chen, Yongqi Huang, Yuning Zhang, and Wanli Ouyang.  Merging vision transformers from different tasks and domains.  _arXiv preprint arXiv:2312.16240_ , 2023.
  * Yu et al. (2023) Le Yu, Bowen Yu, Haiyang Yu, Fei Huang, and Yongbin Li.  Language models are super mario: Absorbing abilities from homologous models as a free lunch.  _arXiv preprint arXiv:2311.03099_ , 2023.
  * Yu et al. (2020) Tianhe Yu, Saurabh Kumar, Abhishek Gupta, Sergey Levine, Karol Hausman, and Chelsea Finn.  Gradient surgery for multi-task learning.  _NeurIPS_ , 33:5824–5836, 2020.
  * Yue et al. (2023) Xiang Yue, Xingwei Qu, Ge Zhang, Yao Fu, Wenhao Huang, Huan Sun, Yu Su, and Wenhu Chen.  Mammoth: Building math generalist models through hybrid instruction tuning.  _arXiv preprint arXiv:2309.05653_ , 2023.
  * Yue et al. (2024) Xiang Yue, Tuney Zheng, Ge Zhang, and Wenhu Chen.  Mammoth2: Scaling instructions from the web.  _arXiv preprint arXiv:2405.03548_ , 2024.
  * Zhang & Yang (2021) Yu Zhang and Qiang Yang.  A survey on multi-task learning.  _IEEE Transactions on Knowledge and Data Engineering_ , 34(12):5586–5609, 2021.
  * Zheng et al. (2024) Chujie Zheng, Ziqi Wang, Heng Ji, Minlie Huang, and Nanyun Peng.  Weak-to-strong extrapolation expedites alignment.  _arXiv preprint arXiv:2404.16792_ , 2024.
  * Zhong et al. (2020) Haoxi Zhong, Chaojun Xiao, Cunchao Tu, Tianyang Zhang, Zhiyuan Liu, and Maosong Sun.  Jec-qa: a legal-domain question answering dataset.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 34, pp. 9701–9708, 2020.
  * Zhong et al. (2023) Wanjun Zhong, Ruixiang Cui, Yiduo Guo, Yaobo Liang, Shuai Lu, Yanlin Wang, Amin Saied, Weizhu Chen, and Nan Duan.  Agieval: A human-centric benchmark for evaluating foundation models.  _arXiv preprint arXiv:2304.06364_ , 2023.

## Appendix

##  Appendix A Proof

###  A.1 Proof of Lemma 4

Using Taylor expansion for $\mathcal{L}(\bm{\theta}_{\textup{final}},\bm{x})$
at $\bm{\theta}_{0}$:

$\displaystyle\mathcal{L}(\bm{\theta}_{\textup{final}},\bm{x})=$
$\displaystyle\mathcal{L}_{t}(\sum_{k=1}^{n}\lambda_{k}(\bm{\theta}_{k}-\bm{\theta}_{0})+\bm{\theta}_{0},\bm{x}_{t})$
(15) $\displaystyle=$
$\displaystyle\mathcal{L}_{t}(\bm{h}_{t}+\bm{\theta}_{t},\bm{x}_{t})$ (16)
$\displaystyle=$
$\displaystyle\mathcal{L}_{t}(\bm{\theta}_{t},\bm{x}_{t})+\nabla\mathcal{L}_{t}(\bm{\theta}_{t},\bm{x}_{t})\bm{h}_{t}+\frac{1}{2}\bm{h}_{t}^{\top}\left(\int_{0}^{1}\nabla^{2}\mathcal{L}_{t}(\gamma_{t}(\beta))d\beta\right)\bm{h}_{t}$

where
$\gamma_{t}(\beta)=\bm{\theta}_{t}+\beta(\bm{\theta}_{\text{final}}-\bm{\theta}_{t})$
and $\bm{h}_{t}$ is the linear combination of $\bm{\lambda}$ and
$\bm{\theta}$:

$\displaystyle\bm{h}_{t}=\sum_{k\neq
t}\lambda_{k}(\bm{\theta}_{k}-\bm{\theta}_{0})-(1-\lambda_{t})(\bm{\theta}_{t}-\bm{\theta}_{0})$
(17)

Because the $\bm{\theta}_{t}$ is fine-tuned using loss $\mathcal{L}_{t}$, the
gradient of $\mathcal{L}_{t}$ at $\bm{\theta}_{t}$ is zero, and the first
order expansion is 0. Substituting Eq. A.1 to Eq. 3, we have:

$\displaystyle\textup{TLD}_{t}$
$\displaystyle=\mathcal{L}_{t}(\bm{\theta}_{\textup{final}},\bm{x}_{t})-\mathcal{L}_{t}(\bm{\theta}_{t},\bm{x}_{t})$
(18)
$\displaystyle=\frac{1}{2}\bm{h}_{t}^{\top}\left(\int_{0}^{1}\nabla^{2}\mathcal{L}_{t}(\gamma_{t}(\beta))d\beta\right)\bm{h}_{t}$
(19)

Thus, we have completed the proof.

###  A.2 Proof of Theorem 7

Before starting the proof, we first introduce a lemma:

######  Lemma 11.

Under the Property. 5, the task vector is linearly with the gradient.

$\displaystyle\delta_{t}(\bm{\theta}_{t}-\bm{\theta}_{0})=\nabla_{\bm{\theta}_{0}}f(\bm{x},\bm{\theta}_{0})$
(20)

Proof: For gradient descent, we have:

$\displaystyle\bm{\theta}_{t}-\bm{\theta}_{0}$
$\displaystyle=\sum_{i=1}^{n}lr_{i}\nabla\mathcal{L}_{t}^{i}=\sum_{i=1}^{n}lr_{i}\frac{\partial\mathcal{L}_{t}^{i}}{\partial
f}\nabla{f}_{i}$ (21)

where $lr_{i}$ and $\nabla\mathcal{L}_{t}^{i}$ and $\nabla f_{i}$ is the
learning rate, gradient loss, gradient of $f$ at step i. From Property 5, we
can see that the fine-tuning process of $f$ occurs in the linear regime, which
indicates that the first order derivative in the task vector direction is an
constant. We derivative at $\bm{\theta}_{t}$:

$\displaystyle\nabla_{\bm{\theta}_{t}}f(\bm{x},\bm{\theta}_{t})=\nabla_{\bm{\theta}_{0}}f(\bm{x},\bm{\theta}_{0})$
(22)

Thus, we substitute all the gradient of $f_{i}$ using
$\nabla_{\bm{\theta}_{0}}f(\bm{x},\bm{\theta}_{0})$:

$\displaystyle\delta_{t}(\bm{\theta}_{t}-\bm{\theta}_{0})=\nabla_{\bm{\theta}_{0}}f(\bm{x},\bm{\theta}_{0})$
(23)
$\text{where}\quad\frac{1}{\delta_{t}}=\sum_{i=1}^{n}lr_{i}\frac{\partial\mathcal{L}_{t}^{i}}{\partial
f}$

Thus, we have completed the proof of the Lemma.

For the of loss function, using Property 5 we have:

$\displaystyle\mathcal{L}_{t}(\bm{\theta}_{t},\bm{x}_{t})=$
$\displaystyle\frac{1}{2}\left\|f(\bm{x}_{t},\bm{\theta}_{t})-y\right\|^{2}=\frac{1}{2}\|(\bm{\theta}_{t}-\bm{\theta}_{0})^{\top}\nabla
f(\bm{x}_{t};\bm{\theta}_{0})+C_{0}\|^{2}$ (24)

For the Hessian of loss function, it can be represented as:

$\displaystyle\nabla^{2}_{\bm{\theta}_{t}}\mathcal{L}_{t}$
$\displaystyle=\nabla_{\bm{\theta}_{0}}f(\bm{x}_{t};\bm{\theta}_{0})\nabla^{\top}_{\bm{\theta}_{0}}f(\bm{x}_{t};\bm{\theta}_{0})$
(25)

Using Eq. 25 the $\textup{TLD}_{t}$ can be represented as:

$\displaystyle 2\textup{TLD}_{t}=$
$\displaystyle\bm{h}_{t}^{\top}\left(\int_{0}^{1}\nabla^{2}\mathcal{L}_{t}(\gamma_{t}(\beta))d\beta\right)\bm{h}_{t}$
(26) $\displaystyle=$
$\displaystyle\bm{h}_{t}^{\top}\left(\nabla^{2}\mathcal{L}_{t}(\tilde{\bm{\theta}})\right)\bm{h}_{t}$
(27) $\displaystyle=$
$\displaystyle\bm{h}_{t}^{\top}\left(\nabla_{\bm{\theta}_{0}}f(\bm{\theta}_{0},\bm{x}_{t})\nabla_{\bm{\theta}_{0}}f^{\top}(\bm{\theta}_{0},\bm{x}_{t})\right)\bm{h}_{t}$
(28) $\displaystyle=$
$\displaystyle\text{tr}\bigg{\\{}\bm{h}_{t}^{\top}\left(\nabla_{\bm{\theta}_{0}}f(\bm{\theta}_{0},\bm{x}_{t})\nabla_{\bm{\theta}_{0}}f^{\top}(\bm{\theta}_{0},\bm{x}_{t})\right)\bm{h}_{t}\bigg{\\}}$
(29) $\displaystyle\leq$
$\displaystyle\text{tr}(\bm{h}\bm{h}^{\top})\text{tr}\left(\nabla_{\bm{\theta}_{0}}f(\bm{\theta}_{0},\bm{x}_{t})\nabla_{\bm{\theta}_{0}}f^{\top}(\bm{\theta}_{0},\bm{x}_{t})\right)$
(30)

For $\text{tr}(\bm{h}\bm{h}^{\top})$, using Property. 6, we have:

$\displaystyle\text{tr}(\bm{h}\bm{h}^{\top})=$
$\displaystyle\left\|\sum_{k\neq
t}^{T}\lambda_{k}(\bm{\theta}_{k}-\bm{\theta}_{0})-(1-\lambda_{t})(\bm{\theta}_{t}-\bm{\theta}_{0})^{\top}\right\|^{2}$
(31) $\displaystyle=$ $\displaystyle\sum_{k\neq t}^{T}\left[\vmathbb{1}_{k\neq
t}(\lambda^{2}_{k})+\vmathbb{1}_{k=t}(1-\lambda^{2}_{k})\right]\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}$
(32) $\displaystyle=$ $\displaystyle\sum_{k\neq
t}^{T}\left[\vmathbb{1}(\lambda^{2}_{k})\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\right]$
(33) $\textup{where}\quad(\lambda^{2}_{k})\vmathbb{1}(k\neq
t)+(1-\lambda^{2}_{k})\vmathbb{1}(k=t):=\vmathbb{1}_{t}(\lambda_{k}^{2}).$

For the second part:
$\text{tr}(\nabla_{\bm{\theta}_{0}}f(\bm{\theta}_{0},\bm{x}_{t})\nabla_{\bm{\theta}_{0}}f^{\top}(\bm{\theta}_{0},\bm{x}_{t}))$,
using Lemma 11 we can have:

$\displaystyle\text{tr}\left(\nabla_{\bm{\theta}_{0}}f(\bm{\theta}_{0},\bm{x}_{t})\nabla_{\bm{\theta}_{0}}f^{\top}(\bm{\theta}_{0},\bm{x}_{t})\right)=\delta_{t}^{2}\|\bm{\theta}_{t}-\bm{\theta}_{0}\|^{2}$
(34)

Thus, for $\textup{TLD}_{t}$ we can upper bound it by

$\displaystyle\textup{TLD}_{t}$
$\displaystyle\leq\frac{\delta_{t}^{2}}{2}\left\|\bm{\theta}_{t}-\bm{\theta}_{0}\right\|_{2}^{2}\bigg{\\{}\sum_{k\neq
t}^{T}\vmathbb{1}_{t}(\lambda_{k}^{2})\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\bigg{\\}}$
(35)

###  A.3 Proof of Theorem 8

By summing Eq.35 from 1 to T, we can complete the proof.

ALD
$\displaystyle\leq\sum_{t=1}^{T}\frac{\delta_{t}^{2}}{2}\left\|\bm{\theta}_{t}-\bm{\theta}_{0}\right\|_{2}^{2}\bigg{\\{}\sum_{k\neq
t}^{T}\vmathbb{1}_{t}(\lambda_{k}^{2})\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\bigg{\\}}$
(36)

###  A.4 Proof of Theorem 9

First, for Eq. 36, we have:

ALD
$\displaystyle\leq\frac{\delta_{t}^{2}}{2}\sum_{t=1}^{T}\left\|\bm{\theta}_{t}-\bm{\theta}_{0}\right\|_{2}^{2}\bigg{\\{}\sum_{k\neq
t}^{T}\vmathbb{1}_{t}(\lambda_{k}^{2})\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\bigg{\\}}$
(37)

where $\delta_{0}=\max\\{\delta_{i}\\}$ For Eq. 37, it is easy to verify that
the terms containing $\lambda_{t}$ can be represented as:

$\displaystyle\textup{ALD}_{\lambda_{t}}=$
$\displaystyle\frac{\delta_{t}^{2}}{2}\|\bm{\theta}_{t}-\bm{\theta}_{0}\|^{2}\left[\sum_{k=1}^{T}\vmathbb{1}_{t}(\lambda)\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\right]$
(38)

Thus, the ALD can be upper bounded by

$\displaystyle\textup{ALD}\leq\sum_{t=1}^{T}\textup{ALD}_{\lambda_{t}}$ (39)

###  A.5 Proof of Theorem 10

Because each $\text{ALD}_{\lambda_{t}}$ does not contain other scaling
coefficients. We can solve each optimal $\lambda_{t}$ from
$\text{ALD}_{\lambda_{t}}$:

$\displaystyle\lambda_{t}$
$\displaystyle=\mathop{\arg\min}_{\lambda_{t}}\frac{\delta_{0}^{2}}{2}\|\bm{\theta}_{t}-\bm{\theta}_{0}\|^{2}\left[\sum_{k=1}^{T}\vmathbb{1}_{t}(\lambda)\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\right]$
(40)
$\displaystyle=\mathop{\arg\min}_{\lambda_{t}}\|\bm{\theta}_{t}-\bm{\theta}_{0}\|^{2}\left[\sum_{k=1}^{T}\vmathbb{1}_{t}(\lambda)\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}\right]$
(41)

The RHS of the above equation is quadratic on $\lambda_{t}$ and and the
optimal solution for $\lambda_{t}$ is:

$\displaystyle\lambda_{t}=\frac{\|\bm{\theta}_{t}-\bm{\theta}_{0}\|^{2}}{\sum_{k=1}^{n}\|\bm{\theta}_{k}-\bm{\theta}_{0}\|^{2}}$
(42)

##  Appendix B Details of Models and Datasets

Table 7 shows the versions and correspondence with pre-trained backbones of
fine-tuned LLMs. Table 8 shows the details of the datasets we use in our
paper.

##  Appendix C Infra and hardware details

We use PyTorch as the deep learning framework. We merge and evaluate the
neural networks using A100 GPUs.

##  Appendix D Hyper-parameter Setting

For both DARE and TIES-Merging, the density of 0.55 is used, and the open-
source tool MergeKit111MergeKit is employed for the merging process.

##  Appendix E Details of different Methods

We give a detailed comparison of the current merging method below from the
perspective of extra data information, time complexity, and optimal
performance. The time complexity for forward and backward processes is denoted
as FW and BP. For RegMean, it requires the inner product data matrices for
layer input to calculate the updated parameters. It only requires a forward
process, but loading all the inner products of the layer input matrix requires
$\mathcal{O}(\theta^{2})$ memory. For Fisher merge, it also requires the data
to calculate the Fisher Matrix, which requires the forward process to
calculate the Fisher matrix and $\mathcal{O}(\theta^{2})$ memory to store the
Fisher matrix. Grid-search Task Arithmetic (G-Task Arithmetic) requires
$\mathcal{O}(G\textsuperscript{T}\times\mathcal{T}_{\text{FW}})$ forward
process to evaluate, where G is the grid number (G = 100 means 100 girds from
0 to 1) and T is the number of tasks. The space complexity is also equal to
the memory requirement of the forward process. For Adamerging, it
simultaneously loads T LLMs to optimize, whose time complexity is
$\mathcal{O}(\mathcal{T}_{\text{BP}})$ and space complexity is:
$\mathcal{O}(\mathcal{S}_{\text{BP}}\times T)$. For weight average, task
arithmetic, and MetaGPT, they all do not need extra data information, which is
model exclusive. Their time and space complexity is $\mathcal{O}(1)$ and
$\mathcal{O}(n)$, but only our MetaGPT achieves optimal performance.

Table 6: Extra data information requirement, time and space complexity, and
optimally of current methods. The time complexity for forward process and back
propagation are denote by $\mathcal{T}_{\text{FW}},\mathcal{T}_{\text{BP}}$.
The space complexity for forward process and back propagation are denote by
$\mathcal{S}_{\text{FW}},\mathcal{S}_{\text{BP}}$. T is the number of task,
$\theta$ is the number of parameters and G is the grid number (G = 100 means
100 girds from 0 to 1).

| Extra Data Info | Time Complexity | Space Complexity | Optimal | Apply to LLMs
---|---|---|---|---|---
RegMean | ✓ | $\mathcal{O}(\mathcal{T}_{\text{FW}})$ | $\mathcal{O}(\theta^{2})$ | ✓ | ✗
Fisher Merge | ✓ | $\mathcal{O}(\mathcal{T}_{\text{FW}})$ | $\mathcal{O}(\theta^{2})$ | ✓ | ✗
G-Task Arithmetic | ✓ | $\mathcal{O}(G\textsuperscript{T}\times\mathcal{T}_{\text{FW}})$ | $\mathcal{O}(\mathcal{S}_{\text{FW}})$ | ✓ | ✗
AdaMerging | ✓ | $\mathcal{O}(\mathcal{T}_{\text{BP}})$ | $\mathcal{O}(\mathcal{S}_{\text{BP}}\times T)$ | ✓ | ✗
Task Arithmetic | ✗ | $\mathcal{O}(1)$ | $\mathcal{O}(\theta)$ | ✗ | ✓
Weight Average | ✗ | $\mathcal{O}(1)$ | $\mathcal{O}(\theta)$ | ✗ | ✓
MetaGPT | ✗ | $\mathcal{O}(1)$ | $\mathcal{O}(\theta)$ | ✓ | ✓

Table 7: Details of datasets we used for our evaluation.

Dataset | Number of Training Examples | Number of Validation Examples | Number of Testing Examples | Evaluate Metric
---|---|---|---|---
WinoGrande | 9,248 | 1,267 | 1,767 | 0-shot accuracy
AGIEval | N/A | N/A | 8,062 | 5-shot accuracy
GSM8k | 7,473 | N/A | 1,319 | 4-shot accuracy
Math | 7,500 | N/A | 1,500 | 4-shot accuracy
MBPP | 374 | 30 | 500 | 3-shot accuracy
HumanEval | N/A | N/A | 164 | 0-shot accuracy
JEC-QA | N/A | N/A | 26,365 | 5-shot accuracy
FinancelQ | N/A | N/A | 7,173 | 5-shot accuracy
MedQA | N/A | N/A | 61,097 | 5-shot accuracy

Table 8: Details of models we used for our evaluation. Pre-trained Model | Task | Fine-tuned-Models
---|---|---
LLaMA-2-7b | General Knowledge | meta-llama/Llama-2-7b-chat-hf
Mathematical Reasoning | TIGER-Lab/MAmmoTH-7B
Code Generating | mrm8488/llama-2-coder-7b
Chinese | hfl/chinese-llama-2-7b
Spanish | clibrain/Llama-2-7b-ft-instruct-es
Japanese | elyza/ELYZA-japanese-Llama-2-7b
Mistral-7b | General Knowledge | mistralai/Mistral-7B-Instruct-v0.2
Mathematical Reasoning | TIGER-Lab/MAmmoTH2-7B
Code Generating | Nondzu/Mistral-7B-codealpaca-lora
LLaMA-2-13b | General Knowledge | meta-llama/Llama-2-13b-chat-hf
Mathematical Reasoning | TIGER-Lab/MAmmoTH-13B
Code Generating | emre/llama-2-13b-code-chat

##  Appendix F Merging Checkpoints of the Pre-trained Model

We further conducted experiments using intermediate checkpoints from the open-
source model Baichuan2-7B Yang et al. (2023a). Surprisingly, we found that
applying our MetaGPTto pre-trained model also lead to performance improvements
across multiple tasks. We tested our approach on MMLU, CMMLU, GSM8k, MATH, and
MBPP evaluations. The results showed that the merged model achieved an average
score of 29.376, which surpasses the best individual checkpoint score of
29.334 and the weight average score of 27.800. Moreover, our approach
consistently outperformed the baseline method of Weight Average. These
findings underscore the robustness and effectiveness of our methodology in
enhancing model performance beyond standard averaging techniques.

Table 9: Performance of different checkpoints.

# of Tokens | MMLU | CMMLU | GSM8k | MATH | MBPP | Avg
---|---|---|---|---|---|---
1700B | 49.95 | 51.17 | 16.15 | 4.0 | 18.0 | 27.854
1710B | 48.67 | 51.61 | 16.53 | 3.8 | 17.8 | 27.682
1720B | 49.85 | 51.37 | 17.06 | 2.8 | 18.8 | 27.976
1730B | 49.91 | 50.24 | 17.44 | 4.2 | 19.2 | 28.198
1740B | 50.99 | 49.64 | 19.64 | 3.8 | 22.6 | 29.334
1750B | 49.50 | 51.20 | 17.13 | 4.6 | 19.2 | 28.326
1760B | 49.02 | 48.52 | 17.29 | 4.2 | 17.8 | 27.366
1770B | 51.48 | 49.74 | 17.36 | 4.4 | 20.4 | 28.676
Weight Average | 48.95 | 50.01 | 17.44 | 4.4 | 18.2 | 27.800
MetaGPT(ours) | 50.61 | 52.69 | 19.18 | 4.4 | 20.0 | 29.376


# Optimal personalized therapies in colon-cancer induced immune response using
a Fokker-Planck framework

Souvik Roy1, Suvra Pal$\,{}^{2}$,

1 Department of Mathematics, The University of Texas at Arlington, Arlington,
TX 76019, USA
2 Department of Mathematics, The University of Texas at Arlington, Arlington,
TX 76019, USA

(Email<EMAIL_ADDRESS>1<EMAIL_ADDRESS>2)

###### Abstract

In this paper, a new stochastic framework to determine optimal combination
therapies in colon cancer-induced immune response is presented. The dynamics
of colon cancer is described through an Itö stochastic process, whose
probability density function evolution is governed by the Fokker-Planck
equation. An open-loop control optimization problem is proposed to determine
the optimal combination therapies. Numerical results with combination
therapies comprising of the chemotherapy drug Doxorubicin and immunotherapy
drug IL-2 validate the proposed framework.

Keywords: Fokker-Planck optimization, non-linear conjugate gradient,
immunotherapy, chemotherapy.

##  3.1 Introduction

Colon cancer is a leading cause of global cancer related deaths [11]. The lack
of early symptoms forces the detection of colon cancer to take place at the
metastatic phase of the cancer [29]. Thus, it becomes important to devise fast
and accurate treatment strategies for cure. In this context, combination
therapies have been clinically shown to be an effective strategy for combating
cancer, in comparison to mono-therapy (see [18]). Conventional mono-
therapeutic techniques are indiscriminate in choosing actively growing cells
that leads to the death of not only cancerous cells but also healthy cells.
E.g., chemotherapy drugs can be toxic and lead to multiple side effects and
risks, by weakening the patient’s immune system targeting the bone-marrow
cells [22]. This leads to increased susceptibility to secondary infections.
However, with combination therapies, even though it might still be toxic if
there is a presence of a chemotherapy drug, the toxicity effect is
significantly diminished due to different targets being affected. Moreover,
since the effect of combination therapies work in synergy, lower dosages of
the individual drugs in the combination therapy are required for treatment,
which further reduces the toxic effects [19].

To develop optimal combination therapies in colon cancer, it is important to
understand relevant biomarkers that govern the progression of the cancer [8].
One such biomarker is the immune response to the cancerous cells. It has been
observed that onsite immune reactions at the tumor locations and prognosis are
highly correlated and is independent of the size of the tumor [12].
Furthermore, a high expression of various immune pathways, like Th1, Th17, are
associated with a poor prognosis or prolonged disease-free survival for
patients with colon cancer [32]. Thus, it becomes important to develop optimal
personalized treatments for colon cancer patients, taking into account the
various types of immune cells, their numbers and interaction with the cancer
cells and each other. Since drugs for treatment of colon cancer, like
chemotherapeutic drugs, immunotherapies, are quite costly, it is expensive to
perform in-vitro and vivo experimental studies for testing effects of such
drugs. An alternative cost effective option is to develop computational
frameworks for testing optimal combination drug dosages [7]. This can be done
using pharmacokinetic models, that are given by a set of ordinary or partial
differential equations (ODE or PDE) .

There are several dynamical models that have been developed to describe the
immune response and interactions in colon cancer and associated therapies. In
[4], a combined compartmental model with that of irinotecan is used to
describe the physiology of colon cancer. A set of ODEs were used to represent
a 3D structure of colon cancer in [14]. In [15], the authors use describe the
phenomenon of tumorigenesis in colon cancer using mathematical modeling. The
authors in [17], describe the initiation of colon cancer and its link with
colitis using a mathematical framework. In [23], the authors use a
pharmacokinetic cellular automata model to incorporate the cytotoxic effects
of chemotherapy drugs. In [13], dynamical systems are used to describe
multiple pathways in colon cancer. For a detailed review of various
mathematical models for colon cancer-induced immune response, we refer the
reader to [3].

We contribute to the field of pharmacokinetic cancer research by presenting an
effective approach to develop personalized therapies for colon cancer-induced
immune response. The starting point of this estimation process is the
dynamical model for colon cancer-induced immune response, given in [10]. The
model describes the evolution of four variables: the tumor cell count,
concentration of the Natural Killer cells, concentration of the CD$8^{+}$
cells, concentration of the remaining lymphocytes. We extend the dynamical
model in [10] to an Itô stochastic ODE system that takes into account the
randomness of the colon cancer-induced immune response dynamics. We also model
a combination therapy comprising of a chemotherapy drug and an immunotherapy
drug as a control in the stochastic ODE system. The goal is to obtain an
optimal combination therapy strategy that drives the cancer-induced immune
response state. Since solving an optimal control problem using stochastic
state equations is difficult, we use a more convenient framework for obtaining
optimal combination therapies through the Fokker-Planck (FP) equations, that
governs the evolution of the joint probability density function (PDF)
associated to the random variables in stochastic process. Such FP control
frameworks have been used for problems arising in control of collective and
crowd motion [24, 25], investigating pedestrian motion from a game theoretic
perspective [26], reconstructing cell-membrane potentials [2], mean field
control problems [5], controlling production of Subtilin [31]. Very recently,
in [28, 27], FP frameworks were used for parameter estimation in colon cancer-
induced dynamical systems. Till date, this is the first work that considers
the FP framework to devise optimal combination treatment strategies for colon
cancer-induced immune response.

In the next section, we describe a a FP control framework for devising optimal
treatments in colon cancer-induced immune response systems. Section 3.3 is
concerned with the theoretical properties of the FP optimization problems. In
Section 3.4, we describe the numerical discretization of the optimality
system. In Section 3.5, we obtain the optimal combination therapies,
comprising of Doxorubicin and IL-2 for two types of simulated colon cancer
patients, and show the correspondence of the results with experimental
findings. We end with a section of conclusions.

##  3.2 The control framework based on Fokker-Planck equations

The interaction between immune cells and colon cancer growth in a patient can
be modeled using a coupled system of ODEs. The starting point is the set of
equations given in the paper by dePillis [10]. The following are the variables
associated with different types of cell populations that appear as model
variables

  1. 1.

$T(\tau)$\- tumor cell population number (cells).

  2. 2.

$N(\tau)$\- natural killer (NK) cells concentration per liter of blood
(cells/L).

  3. 3.

$L(\tau)$\- cytotoxic T lymphocytes (CD8+) concentration per liter of blood
(cells/L).

  4. 4.

$C(\tau)$\- other lymphocytes concentration per liter of blood (cells/L).

The governing system of ODEs, representing the dynamics of the above defined
cell populations are given as follows

$\displaystyle\dfrac{dT}{d\tau}=aT-abT^{2}-DT-
cNT-\alpha_{1}u_{1}(t)T,~{}T(0)=T_{0},$ (3.1)
$\displaystyle\dfrac{dN}{d\tau}=eC-fN-
pNT-\alpha_{2}u_{1}(t)N+\beta_{2}u_{2}(t)N,~{}N(0)=N_{0},$
$\displaystyle\dfrac{dL}{d\tau}=mL+j\dfrac{T}{k+T}L-qLT+(r_{1}N+r_{2}C)T-\alpha_{3}u_{1}(t)L+\beta_{3}u_{2}(t)L,~{}L(0)=L_{0},$
$\displaystyle\dfrac{dC}{d\tau}=\alpha-\beta
C-\alpha_{4}u_{1}(t)C,~{}C(0)=C_{0},$

where $D=d\dfrac{1}{s(T/L)^{l}+1}$, $T_{0},N_{0},L_{0},C_{0}$ represent the
initial conditions for $T,N,L,C$, respectively, and $u_{1},u_{2}$ represent
dosages of Doxorubicin and IL-2, respectively.

The parameters of the system (3.1) are defined in [10]. To stabilize the
solutions of the numerical methods, we non-dimensionalize the ODE system (3.1)
using the following non-dimensionalized state variables and parameters

$\displaystyle\bar{T}=k_{1}T,~{}\bar{N}=k_{2}N,~{}\bar{L}=k_{3}L,~{}\bar{C}=k_{4}C,~{}t=k_{5}\tau,$
(3.2)
$\displaystyle\bar{a}=\dfrac{a}{k_{5}},~{}\bar{b}=\dfrac{b}{k_{1}},~{}\bar{c}=\dfrac{c}{k_{2}k_{5}},~{}\bar{d}=\dfrac{d}{k_{5}},~{}\bar{e}=\dfrac{ek_{2}}{k_{4}k_{5}},~{}\bar{f}=\dfrac{f}{k_{5}},~{}\bar{p}=\dfrac{10^{10}p}{k_{1}k_{5}},~{}\bar{m}=\dfrac{m}{k_{5}},~{}\bar{j}=\dfrac{j}{k_{5}},~{}\bar{k}=k_{1}k,$
$\displaystyle\bar{q}=\dfrac{10^{8}q}{k_{5}k_{1}},~{}\bar{r}_{1}=\dfrac{r_{1}k_{3}}{k_{1}k_{2}k_{5}},~{}\bar{r}_{2}=\dfrac{r_{2}k_{3}}{k_{1}k_{4}k_{5}},~{}\bar{s}=250s,~{}\bar{\alpha}=\dfrac{\alpha
k_{4}}{k_{5}},~{}\bar{\beta}=\dfrac{\beta}{k_{5}}~{}\bar{l}=l,$
$\displaystyle\bar{\alpha}_{i}=\dfrac{\alpha_{i}}{k_{5}},~{}i=1,\cdots,4,$
$\displaystyle\bar{\beta}_{i}=\dfrac{\beta_{i}}{k_{5}}~{}i=1,2.$

Then the transformed non-dimensional ODE system is given as follows

$\displaystyle\dfrac{d\bar{T}}{dt}=\bar{a}\bar{T}-\bar{a}\bar{b}\bar{T}^{2}-\bar{D}\bar{T}-\bar{c}\bar{N}\bar{T}-\bar{\alpha}_{1}\bar{u}_{1}(t)\bar{T},~{}\bar{T}(0)=\bar{T}_{0}$
(3.3)
$\displaystyle\dfrac{d\bar{N}}{dt}=\bar{e}\bar{C}-\bar{f}\bar{N}-10^{-10}\bar{p}\bar{N}\bar{T}-\bar{\alpha}_{2}\bar{u}_{1}(t)\bar{N}+\bar{\beta}_{2}\bar{u}_{2}(t)\bar{N},~{}\bar{N}(0)=\bar{N}_{0}$
$\displaystyle\dfrac{d\bar{L}}{dt}=\bar{m}\bar{L}+\bar{j}\dfrac{\bar{T}}{\bar{k}+\bar{T}}\bar{L}-10^{-8}\bar{q}\bar{L}\bar{T}+(\bar{r}_{1}\bar{N}+\bar{r}_{2}\bar{C})\bar{T}-\bar{\alpha}_{3}\bar{u}_{1}(t)\bar{L}+\bar{\beta}_{3}\bar{u}_{2}(t)\bar{L},~{}\bar{L}(0)=\bar{L}_{0}$
$\displaystyle\dfrac{d\bar{C}}{dt}=\bar{\alpha}-\bar{\beta}\bar{C}-\bar{\alpha}_{4}\bar{u}_{1}(t)\bar{C},~{}\bar{C}(0)=\bar{C}_{0},$

where $\bar{D}=\bar{d}\dfrac{(\bar{L}/\bar{T})^{\bar{l}}}{4\bar{s}\cdot
10^{-3}\cdot(k_{1}/k_{3})^{\bar{l}}+(\bar{L}/\bar{T})^{\bar{l}}}$.

The compact form of the aforementioned system of ODEs, given in (3.3), is as
follows:

$\displaystyle\dfrac{d\bm{X}}{dt}=\bm{F}(\bm{X},\bm{U}),$ (3.4)
$\displaystyle\bm{X}(0)=\bm{X}_{0},$

where $\bm{X}(t)=(\bar{T}(t),\bar{N}(t),\bar{L}(t),\bar{C}(t))^{T}$.

We extend the ODE system (3.3) to include stochasticity present in the
dynamics. For this purpose, we consider the Itô stochastic differential
equation corresponding to (3.3)

$\displaystyle
d\bar{T}=(\bar{a}\bar{T}(1-\bar{b}\bar{T})-\bar{c}\bar{N}\bar{T}-\bar{D}\bar{T}-\bar{\alpha}_{1}\bar{u}_{1}(t)\bar{T})~{}dt+\sigma_{1}(\bar{T})~{}dW_{1}(t),~{}\bar{T}(0)=\bar{T}_{0},$
(3.5) $\displaystyle
d\bar{N}=(\bar{e}\bar{C}-\bar{f}\bar{N}-10^{-10}\bar{p}\bar{N}\bar{T}-\bar{\alpha}_{2}\bar{u}_{1}(t)\bar{N}+\bar{\beta}_{2}\bar{u}_{2}(t)\bar{N})~{}dt+\sigma_{2}(\bar{N})~{}dW_{2}(t),~{}\bar{N}(0)=\bar{N}_{0},$
$\displaystyle
d\bar{L}=(\bar{m}\bar{L}+\bar{j}\dfrac{\bar{T}}{\bar{k}+\bar{T}}\bar{L}-10^{-8}\bar{q}\bar{L}\bar{T}+(\bar{r}_{1}\bar{N}+\bar{r}_{2}\bar{C})\bar{T}-\bar{\alpha}_{3}\bar{u}_{1}(t)\bar{L}+\bar{\beta}_{3}\bar{u}_{2}(t)\bar{L})~{}dt+\sigma_{3}(\bar{L})~{}dW_{3}(t),$
$\displaystyle\bar{L}(0)=\bar{L}_{0},$ $\displaystyle
d\bar{C}=(\bar{\alpha}-\bar{\beta}\bar{C}-\bar{\alpha}_{4}\bar{u}_{1}(t)\bar{C})~{}dt+\sigma_{4}(\bar{C})~{}dW_{4}(t),~{}\bar{C}(0)=\bar{C}_{0},$

where $dW_{i},~{}i=1,2,3,4$ are one-dimensional Wiener processes and
$\sigma_{i},~{}i=1,2,3,4$ are positive constants. The equation (3.5) can be
written using a compact notation as follows

$\displaystyle
d\bm{X}=\bm{F}(\bm{X},\bm{U})~{}dt+\bm{\sigma}(\bm{X})~{}d\bm{W}(t),$ (3.6)
$\displaystyle\bm{X}(0)=\bm{X}_{0},$

where

$d\bm{W}(t)=\begin{pmatrix}dW_{1}(t)&dW_{2}(t)&dW_{3}(t)&dW_{4}(t)\\\
\end{pmatrix}^{T}$

is a four-dimensional Wiener process vector with stochastically independent
components and

$\bm{\sigma}=diag\begin{pmatrix}\sigma_{1}&\sigma_{2}&\sigma_{3}&\sigma_{4}\end{pmatrix}$

is the dispersion matrix.

We now describe the PDF of the stochastic process (3.6), confined in a
Lipschitz domain $\Omega$, by virtue of a reflecting barrier on
$\partial\Omega$. This is motivated by the maximum cell carrying capacity.
Thus,
$X(t)\in\Omega\subset{\mathbb{R}}^{4}_{+}=\\{x\in\mathbb{R}^{4}:x_{i}\geq
0,~{}i=1,2,3,4\\}$, Let $x=(x_{1},x_{2},x_{3},x_{4})^{T}$. Define $f(x,t)$ as
the PDF for the stochastic process described by (3.6), i.e., $f(x,t)$ is the
probability of $\bm{X}(t)$ assuming the value $x$ at time $t$. Then the PDF of
$\bm{X}(t)$ evolves through the following Fokker-Planck (FP) equations

$\displaystyle\partial_{t}f(x,t)+\nabla\cdot(\bm{F}(x,\bm{U})~{}f(x,t))=\frac{1}{2}\nabla\cdot(\bm{\sigma}^{2}(x)\nabla
f(x,t)),$ (3.7) $\displaystyle f(x,0)=f_{0}(x),$

where $f_{0}(x)$ is non-negative with mass equals one, and
$\bm{U}(t)=(\bar{u}_{1}(t),\bar{u}_{2}(t))$ in the admissible set

$U_{ad}=\\{\bm{U}\in L^{2}([0,T]):0\leq\bar{u}_{i}(t)\leq
D_{i},~{}D_{i}>0,~{}\forall t\in[0,T],~{}i=1,2\\}$

Here, the FP domain is $Q=\Omega\times(0,T_{f})$, where $T_{f}$ is the final
time and $f_{0}(x)$ represents the distribution of the initial state $X_{0}$
of the process. The FP equation (3.7) is associated with the no-flux boundary
conditions. To describe this, we write (3.7) in flux form as

$\partial_{t}f(x,t)-\nabla\cdot\mathcal{H}=0,\qquad f(x,0)=f_{0}(x),$ (3.8)

where the components of the flux $\mathcal{H}$ is given as follows

$\mathcal{H}_{j}(x,t;f)=\frac{\sigma_{j}^{2}}{2}\partial_{x_{j}}f-\bm{F}_{j}(x,\bm{U})f,~{}j=1,2,3,4.$
(3.9)

Then the no-flux boundary conditions are

$\mathcal{H}\cdot\hat{n}=0\qquad\mbox{ on }\partial\Omega\times(0,T_{f}),$
(3.10)

with $\hat{n}$ as the unit outward normal on $\partial\Omega$.

To obtain the optimal combination therapy function vector $\bm{U}$, we solve
the following optimization problem

$\displaystyle\bm{U}^{*}=\operatorname*{arg\,min}_{\bm{U}\in
U_{ad}}J(f,\bm{U}):=$
$\displaystyle\dfrac{\alpha}{2}\int_{Q}(f(x,t)-f^{*}(x,t))^{2}~{}dx+\dfrac{\nu_{1}}{2}\int_{0}^{T}\|u_{1}(t)\|^{2}~{}dt+\dfrac{\nu_{2}}{2}\int_{0}^{T}\|u_{2}(t)\|^{2}~{}dt,$
(3.11)

subject to the FP system (3.7),(3.10), where the desired PDF is $f^{*}(x,t)$.

##  3.3 Theoretical results

In this section, we describe some theoretical results related to the
minimization problem (3.11). One can also find similar results in [1, 24, 25].
For this purpose, we denote the FP system (3.7),(3.10) as
${\mathcal{E}}(f_{0},\bm{U})=0$. The existence and uniqueness of solutions of
(3.7) is given in the following proposition

######  Proposition 3.1.

Assume $f_{0}\in H^{1}(\Omega)$ with $f_{0}$ non-negative, and $\bm{U}\in
U_{ad}$. Then, there exists an unique non-negative solution of
${\mathcal{E}}(f_{0},\bm{U})=0$ given by $f\in
L^{2}([0,T_{f}];H^{1}(\Omega))\cap C([0,T_{f}];L^{2}(\Omega))$.

Under the assumptions of higher regularity of $\partial\Omega$, the boundary
of $\Omega$, one can also obtain $H^{2}(\Omega)$ regularity of the solution of
(3.7) (see [30]). Next we state the conservativeness property of (3.7), which
can be proved using straightforward applications of weak formulation,
integration by parts, and divergence theorem for the flux $\mathcal{H}$.

######  Proposition 3.2.

The FP system given in (3.7),(3.10) is conservative.

The next proposition states and proves the $L^{2}$ stability property of
(3.7).

######  Proposition 3.3.

The FP system (3.7),(3.10) solution, given by $f$, satisfies the following
$L^{2}$ stability property

$\|f(t)\|_{L^{2}(\Omega)}\leq\|f_{0}\|_{L^{2}(\Omega)}\exp\left(\|\bm{\sigma}^{-1}\|^{2}_{2}N^{2}t\right),$
(3.12)

where $N=\sup_{\Omega\times U}|F(x,\bm{\theta})|$.

###### Proof.

Multiplying (3.7) with the test function $\psi=f(\cdot,t)$ and an integration
by parts gives

$\frac{\partial}{\partial t}\|f(t)\|^{2}_{L^{2}(\Omega)}=-\|\bm{\sigma}\nabla
f(t)\|^{2}_{L^{2}(\Omega)}+2\int\limits_{\Omega}(\bm{F}f(t))\cdot\bm{\sigma}^{-1}\bm{\sigma}\nabla
f(t)\,dx.$ (3.13)

The last term in (3.13) can be estimated using the Young’s inequality,
$2bd\leq kb^{2}+d^{2}/k$ with $k=\|\bm{\sigma}^{-1}\|_{2}$, which is the
$L^{2}$ matrix norm of $\bm{\sigma}^{-1}$. We then obtain the following

$\frac{\partial}{\partial
t}\|f(t)\|^{2}_{L^{2}(\Omega)}\leq\|\bm{\sigma}^{-1}\|^{2}_{2}N^{2}\|f(t)\|^{2}_{L^{2}(\Omega)}.$

Applying the Gronwall’s inequality gives the desired result. ∎

The above results imply that the map $\Lambda:U_{ad}\to
C([0,T_{f}];H^{1}(\Omega))$ given by $f=\Lambda(\bm{U})$, is continuous and
Fréchet differentiable. The main theoretical result in this work, which shows
that there exists an optimal open-loop control $\bm{U}^{*}$, is given below:

######  Theorem 3.1.

Let $f_{0}\in H^{1}(\Omega)$ and let $J$ be given as in (3.11). Then there
exists $(f^{*},\bm{U}^{*})\in C([0,T_{f}];$ $H^{1}(\Omega))\times U_{ad}$ with
$f^{*}$ being a solution to ${\mathcal{E}}(f_{0},\bm{U}^{*})=0$ and
$\bm{U}^{*}$ minimizing $J$ in $U_{ad}$.

###### Proof.

Since $J$ is bounded, a minimizing sequence $(\bm{U}^{m})$ exists in $U_{ad}$.
Moreover, $J$ is being coercive and sequentially weakly lower semi-continuous
in $U_{ad}$, implies the boundedness of this sequence. Due to the fact that
$U_{ad}$ is a closed and convex subset of a Hilbert space, the sequence
$(\bm{U}^{m})$ contains a convergent subsequence $(\bm{U}^{m_{l}})$ in
$U_{ad}$, such that $\bm{U}^{m_{l}}\rightarrow\bm{U}^{*}$. Correspondingly,
the sequences $(f^{m_{l}}),(\partial_{t}f^{m_{l}})$, where
$f^{m_{l}}=\Lambda(\bm{\theta}^{m_{l}})$, are bounded in
$L^{2}([0,T_{f}];H^{1}(\Omega)),~{}L^{2}([0,T_{f}];H^{-1}(\Omega))$,
respectively. This implies the weak convergence of the sequences to $f^{*}$
and $\partial_{t}f^{*}$, respectively. We next use the compactness result of
Aubin-Lions [16] to obtain strongly convergent subsequence $(f^{m_{k}})$ in
$L^{2}([0,T_{f}],L^{2}(\Omega))$. Thus, the sequence
$(\bm{F}(\bm{U}^{m_{k}})f^{m_{k}})$ in $L^{2}([0,T_{f}],L^{2}(\Omega))$ is
weakly convergent. This implies that $f^{*}=\Lambda(\bm{U}^{*})$, and
$(f^{*},\bm{U}^{*})$ is a minimizer of $J$. ∎

For the minimization problem (3.11), the optimality system can now be written
as

$\displaystyle\partial_{t}f(x,t)+\nabla\cdot(\bm{F}(x,\bm{U})~{}f(x,t))=\frac{1}{2}\nabla\cdot(\bm{\sigma}^{2}\nabla
f(x,t)),\qquad\mbox{ in }\Omega\times(0,T_{f}),$ (FOR) $\displaystyle
f(x,0)=f_{0}(x),\qquad\mbox{ in }\Omega,$
$\displaystyle\mathcal{H}\cdot\hat{n}=0,\qquad\mbox{ on
}\partial\Omega\times(0,T_{f}).$

$\displaystyle-\partial_{t}p(x,t)-f(x,t)($
$\displaystyle\bm{F}(x,\bm{\theta})\cdot\nabla
p(x,t))-\frac{1}{2}\nabla\cdot(\bm{\sigma}^{2}\nabla
p(x,t))=-\alpha(f(x,t)-f_{i}^{*}(x,t)),~{}\mbox{ in }\Omega\times(0,T_{f}),$
(ADJ) $\displaystyle p(x,T_{f})=0,\qquad\mbox{ in }\Omega,$
$\displaystyle\frac{\partial{p}}{\partial{n}}=0,\qquad\mbox{ on
}\partial\Omega\times(0,T_{f}).$

$\Big{\langle}{\beta\bm{U}-\int_{\Omega}\nabla_{\bm{U}}\bm{F}\cdot\nabla
p,\bm{\psi}-\bm{U}}\Big{\rangle}_{L^{2}([0,T])}\geq
0,\qquad\forall\bm{\psi}\in U_{ad}.$ (OPT)

The optimality system comprises of three sets of equations: the forward or
state equation that governs the FP dynamics (FOR), the adjoint equation (ADJ),
and the optimality condition (OPT). In the next section, we describe numerical
schemes to implement the optimality system.

##  3.4 Numerical schemes

We consider the mesh $\\{\Omega_{h}\\}_{h>0}$ given by

$\Omega_{h}=\\{(x_{1},x_{2},x_{3},x_{4})\in\mathbb{R}^{4}:(x_{1i},x_{2j},x_{3k},x_{4l})=(x_{10}+ih,x_{20}+jh,x_{30}+kh,x_{40}+lh)\\},$

where
$(i,j,k,l)\in\\{{0,\ldots,N_{x_{1}}}\\}\times\\{{0,\ldots,N_{x_{2}}}\\}\times\\{{0,\ldots,N_{x_{3}}}\\}\times\\{{0,\ldots,N_{x_{4}}}\\}\cap\Omega$,
and $N_{x_{i}}$ is the number of discretization points along the $i^{th}$
coordinate direction. Define $\delta{t}=T_{f}/N_{t}$ to be the temporal
discretization step, where $N_{t}$ denotes the maximum number of temporal
steps. This gives us the discretized domain for $\Omega$ as follows

$Q_{h,\delta{t}}=\\{{(x_{1i},x_{2j},x_{3k},x_{4l},t_{m}):\,(x_{1i},x_{2j},x_{3k},x_{4l})\in\Omega_{h},~{}t_{m}=m\delta{t},~{}0\leq
m\leq N_{t}}\\}.$

The value of $f(x,t)$ on $Q_{h,\delta{t}}$ is denoted as $f_{i,j}^{m}$. To
solve the forward Fokker-Planck equation (3.7), we use the scheme described in
[27], that is comprised of the spatial discretization using the Chang-Cooper
(CC) method [6]. The temporal derivative discretization is done using the
four-step, alternate direction implicit Douglas-Gunn (DG4) method. The
combined scheme is referred to as the DG4-CC scheme. It has been shown in [27]
that this scheme is positive, conservative, stable, and second order
convergent in the $L^{1}$ norm. Furthermore, we use the temporal D-G scheme
for the time discretization in the first term, one sided finite difference
discretization for the second term, and central difference for the third term
on the left hand side of the adjoint equation (ADJ).

To solve the optimization problem (3.11), we use a projected non-linear
conjugate gradient scheme (PNCG), as described in [28, 27]. Such a scheme has
also been used for solving optimization problems related to controlling
stochastic crowd motion [24, 25], studying avoidance behavior of pedestrians
using game theory [26], and cure rate models [20, 21]. The PNCG scheme is
described below:

######  Algorithm 3.1 (PNCG scheme).

  1. 1.

Input: Initial guess $u_{0}$. Compute $d_{0}=-\nabla\hat{J}(u_{0})_{H^{1}}$.
Set $k=0$ and maximum number of iterations $k_{max}$, tolerance $TOL$.

  2. 2.

If $(k<k_{max})$ do

  3. 3.

Compute $u_{k+1}=P_{U}\left[u_{k}+\alpha_{k}\,d_{k}\right]$, with $\alpha_{k}$
obtained by the Armijo line-search method.

  4. 4.

Evaluate $g_{k+1}=\nabla\hat{J}(u_{k+1})_{H^{1}}$.

  5. 5.

Evaluate
$\beta_{k}^{HG}=\frac{1}{d_{k}^{T}y_{k}}\left({y_{k}-2d_{k}\frac{\|y_{k}\|^{2}}{d_{k}^{T}y_{k}}}\right)^{T}g_{k+1}.$

  6. 6.

Compute $d_{k+1}=-g_{k+1}+\beta_{k}^{HG}d_{k}$.

  7. 7.

Check for $\|u_{k+1}-u_{k}\|_{2}<TOL$. If yes, terminate.

  8. 8.

Update $k=k+1$.

  9. 9.

End if.

##  3.5 Results

This section describes the numerical results for obtaining optimal treatment
strategies with the aforementioned FP framework. We choose our domain
$\Omega=(0,6)^{4}$ and discretize it using $N_{x_{i}}=51$ points for
$i=1,2,3,4$. The final time $T_{f}$ is chosen to be 10 and the maximum number
of time steps $N_{t}$ is chosen to be 200. The values of the constants used in
converting the ODE system (3.1) to its non-dimensional form given in (3.3) are
given as $k_{1}=10^{-10},~{}k_{2}=10^{-5},~{}k_{3}=10^{-7},~{}k_{4}=10^{-8}$,
and $k_{5}=1$. To obtain the target PDF $f^{*}(x,t)$, we simulate the ODE
system (3.3) with $\bar{u}_{1},\bar{u}_{2}$ set to 0 and with the value of the
non-dimensional parameters $(\bar{d},\bar{l},\bar{s})=(2.1,1.1,1.25)$ that
represents a patient with strong immune response [10]. The values of the other
parameters are taken from [10]. After we obtain the trajectories of
$\bar{T},\bar{N},\bar{L},\bar{C}$, we then choose 20 time points $t_{i}$, and
at each point $(\bar{T}(t_{i}),\bar{N}(t_{i}),\bar{L}(t_{i}),\bar{C}(t_{i}))$
assign a Gaussian PDF given by a normal distribution with variance 0.05. We
finally perform a 5D interpolation to obtain the desired PDF $f^{*}(x,t)$.
Based on a statistical analysis of the dataset given in [10], we choose

$\sigma_{i}(x)=0.5(x_{i}^{1.2}+0.001),~{}i=1,2,3,4.$

The regularization parameters are chosen to be
$\alpha=1,~{}\nu_{1}=\nu_{2}=0.01$. The maximum tolerable dosage for
Doxorubicin is taken to be $7$ mg/day and for IL-2 is taken to be $7.2\times
10^{5}$ IU IL-2/l/day, which implies, $D_{1}=7$ and $D_{2}=0.072$.

In Test Case 1, we simulate a tumor patient with a value of
$(T(0),N(0),L(0),C(0))=(2\times 10^{9},10^{5},10^{7},10^{8})$. This implies
$(\bar{T}(0),\bar{N}(0),\bar{L}(0),\bar{C}(0))=(0.2,1,1,1))$. We also choose
the values of the non-dimensional parameters
$(\bar{d},\bar{l},\bar{s})=(1.3,2,10)$ that represents patient with weak
immune system. The goal is to determine optimal drug dosages
$\bar{u}_{i},~{}i=1,2$ such that the tumor profile is given by the target PDF
$f^{*}$.

(a) $T$ (without treatment) (b) $T$ (with treatment) (c) Doxorubicin profile
(d) IL-2 profile

Figure 3.1: Test Case 1: Plots of the mean tumor profiles and drug dosages
over the 10 day period

Figure 3.1 presents the simulation results obtained with our framework. We
observe that without the combination therapy, the mean tumor cell count will
keep on increasing till it reaches the cell carrying capacity, ultimately
leading to patient death. With the combination therapy, the mean tumor cell
count is brought back to diminishing levels. We also note the optimal dosage
patterns for both Doxorubicin and IL-2 over time. Traditionally, Doxorubicin
is administered every 21 days with a dosage of 142.5 mg [9], that translates
to a total dosage of 70 mg for 10 days. However, from Figure 3.1, we note that
the total dosage of Doxorubicin over 10 days is far less than 70 mg. Moreover,
the observed daily dosage of IL-2 from Figure 3.1 is less than the standard
dosage of IL-2 over a day, which is $2.1\times 10^{6}$ IU IL-2/l. This
suggests that optimal combination therapies administered daily leads to lower
total dosages and, thus, lower toxicity effects.

In Test Case 2, we choose the same initial values of
$\bar{T},\bar{N},\bar{L},\bar{C}$, along with the values of the non-
dimensional parameters $(\bar{d},\bar{l},\bar{s})=(1.6,1.4,2)$ that represents
patient with moderately strong immune system.

(a) $T$ (without treatment) (b) $T$ (with treatment) (c) Doxorubicin (d) IL-2

Figure 3.2: Test Case 2: Plots of the mean tumor profiles and drug dosages
over the 10 day period

Figure 3.2 shows the plots of the tumor profiles without and with treatment,
and the optimal daily dosages of Doxorubicin and IL-2. We observe a similar
behavior as Test Case 1. In addition, we note that the dosages are smaller
compared to the patient in Test Case 1, since the immune system is moderately
strong. This shows the effectiveness of our framework in obtaining optimal
dosages in colon cancer.

##  3.6 Conclusion

In this paper, we presented a new stochastic framework to determine optimal
combination therapies in colon cancer-induced immune response. We considered
the tumor and immune system response dynamics as proposed in [10] and extend
it to a stochastic process to account for incorporating randomness in the
dynamics. We characterized the state of the stochastic process using the PDF,
whose evolution is governed by the FP equation. We then solved an optimal
control problem with open loop controls, to obtain the optimal combination
dosages involving chemotherapy and immunotherapy. Numerical results
demonstrate the feasibility of our proposed framework to obtain small dosages
of the combination drugs leading to lower toxicity whilst preserving the
effectiveness to eliminate the tumor.

## Acknowledgments

The authors were supported by the National Institutes of Health (Grant Number:
R21CA242933). S. Roy was also partially supported by the Interdisciplinary
Research Program, University of Texas at Arlington, Grant number: 2021-772.
The content is solely the responsibility of the authors and does not
necessarily represent the official views of the National Institutes of Health.

## References

  * [1] M. Annunziato and A. Borzì.  A fokker–planck control framework for multidimensional stochastic processes.  Journal of Computational and Applied Mathematics, 237(1):487–507, 2013.
  * [2] M. Annunziato and A. Borzì.  A fokker–planck approach to the reconstruction of a cell membrane potential.  SIAM Journal on Scientific Computing, 43(3):B623–B649, 2021.
  * [3] A. Ballesta and J. Clairambault.  Physiologically based mathematical models to optimize therapies against metastatic colorectal cancer: a mini-review.  Current pharmaceutical design, 20(1):37–48, 2014.
  * [4] A. Ballesta, J. Clairambault, S. Dulong, and F. Levi.  A systems biomedicine approach for chronotherapeutics optimization: focus on the anticancer drug irinotecan.  In New Challenges for Cancer Systems Biomedicine, pages 301–327. Springer, 2012.
  * [5] A. Borzì and L. Grüne.  Towards a solution of mean-field control problems using model predictive control.  IFAC-PapersOnLine, 53(2):4973–4978, 2020.
  * [6] J. Chang and G. Cooper.  A practical difference scheme for fokker-planck equations.  Journal of Computational Physics, 6(1):1–16, 1970.
  * [7] Y. Chang, M. Funk, S. Roy, E. Stephenson, S. Choi, H. V. Kojouharov, B. Chen, and Z. Pan.  Developing a mathematical model of intracellular calcium dynamics for evaluating combined anticancer effects of afatinib and rp4010 in esophageal cancer.  International Journal of Molecular Sciences, 23(3):1763, 2022.
  * [8] Y. Chang, S. Roy, and Z. Pan.  Store-operated calcium channels as drug target in gastroesophageal cancers.  Frontiers in pharmacology, 12:944, 2021.
  * [9] L. de Pillis, K. Renee Fister, W. Gu, C. Collins, M. Daub, D. Gross, J. Moore, and B. Preskill.  Mathematical model creation for cancer chemo-immunotherapy.  Computational and Mathematical Methods in Medicine, 10(3):165–184, 2009.
  * [10] L. DePillis, H. Savage, and A. Radunskaya.  Mathematical model of colorectal cancer with monoclonal antibody treatments.  arXiv preprint arXiv:1312.3023, 2013.
  * [11] C. Fitzmaurice, D. Dicker, A. Pain, H. Hamavid, M. Moradi-Lakeh, M. F. MacIntyre, C. Allen, G. Hansen, R. Woodbrook, C. Wolfe, et al.  The global burden of cancer 2013.  JAMA oncology, 1(4):505–527, 2015.
  * [12] J. Galon, A. Costes, F. Sanchez-Cabo, A. Kirilovsky, B. Mlecnik, C. Lagorce-Pagès, M. Tosolini, M. Camus, A. Berger, P. Wind, et al.  Type, density, and location of immune cells within human colorectal tumors predict clinical outcome.  Science, 313(5795):1960–1964, 2006.
  * [13] S. Haupt, A. Zeilmann, A. Ahadova, H. Bläker, M. von Knebel Doeberitz, M. Kloor, and V. Heuveline.  Mathematical modeling of multiple pathways in colorectal carcinogenesis using dynamical systems with kronecker structure.  PLoS computational biology, 17(5):e1008970, 2021.
  * [14] M. D. Johnston, C. M. Edwards, W. F. Bodmer, P. K. Maini, and S. J. Chapman.  Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer.  Proceedings of the National Academy of Sciences, 104(10):4008–4013, 2007.
  * [15] N. L. Komarova, C. Lengauer, B. Vogelstein, and M. A. Nowak.  Dynamics of genetic instability in sporadic and familial colorectal cancer.  Cancer biology & therapy, 1(6):685–692, 2002.
  * [16] J.-L. Lions.  Quelques méthodes de résolution de problemes aux limites non linéaires.  Paris, Dunod-Gauth. Vill., 1969.
  * [17] W.-C. Lo, E. W. Martin Jr, C. L. Hitchcock, and A. Friedman.  Mathematical model of colitis-associated colon cancer.  Journal of theoretical biology, 317:20–29, 2013.
  * [18] R. B. Mokhtari, T. S. Homayouni, N. Baluch, E. Morgatskaya, S. Kumar, B. Das, and H. Yeger.  Combination therapy in combating cancer.  Oncotarget, 8(23):38022, 2017.
  * [19] R. B. Mokhtari, S. Kumar, S. S. Islam, M. Yazdanpanah, K. Adeli, E. Cutz, and H. Yeger.  Combination of carbonic anhydrase inhibitor, acetazolamide, and sulforaphane, reduces the viability and growth of bronchial carcinoid cell lines.  BMC cancer, 13(1):1–18, 2013.
  * [20] S. Pal and S. Roy.  A new non-linear conjugate gradient algorithm for destructive cure rate model and a simulation study: illustration with negative binomial competing risks.  Communications in Statistics-Simulation and Computation, pages 1–15, 2020.
  * [21] S. Pal and S. Roy.  On the estimation of destructive cure rate model: A new study with exponentially weighted poisson competing risks.  Statistica Neerlandica, 75(3):324–342, 2021.
  * [22] A. H. Partridge, H. J. Burstein, and E. P. Winer.  Side effects of chemotherapy and combined chemohormonal therapy in women with early-stage breast cancer.  JNCI Monographs, 2001(30):135–142, 2001.
  * [23] G. G. Powathil, D. J. Adamson, and M. A. Chaplain.  Towards predicting the response of a solid tumour to chemotherapy and radiotherapy treatments: clinical insights from a computational model.  PLoS computational biology, 9(7):e1003120, 2013.
  * [24] S. Roy, M. Annunziato, and A. Borzì.  A fokker–planck feedback control-constrained approach for modelling crowd motion.  Journal of Computational and Theoretical Transport, 45(6):442–458, 2016.
  * [25] S. Roy, M. Annunziato, A. Borzì, and C. Klingenberg.  A fokker–planck approach to control collective motion.  Computational Optimization and Applications, 69(2):423–459, 2018\.
  * [26] S. Roy, A. Borzı, and A. Habbal.  Pedestrian motion constrained by fp-constrained nash games.  Royal Society Open Science, 4:170648, 2017.
  * [27] S. Roy, S. Pal, A. Manoj, S. Kakarla, J. Villegas, and M. Alajmi.  A fokker-planck framework for parameter estimation and sensitivity analysis in colon cancer.  AIP Conference Proceedings (accepted, in press), 2021.
  * [28] S. Roy, Z. Pan, and S. Pal.  A fokker-planck feedback control framework for optimal personalized therapies in colon cancer-induced angiogenesis.  Journal of Mathematical Biology, 84:23, 2022.
  * [29] H. Schmoll, E. Van Cutsem, A. Stein, V. Valentini, B. Glimelius, K. Haustermans, B. Nordlinger, C. Van de Velde, J. Balmana, J. Regula, et al.  Esmo consensus guidelines for management of patients with colon and rectal cancer. a personalized approach to clinical decision making.  Annals of oncology, 23(10):2479–2516, 2012.
  * [30] T. Tao.  Nonlinear dispersive equations: local and global analysis.  Number 106. American Mathematical Soc., 2006.
  * [31] V. Thalhofer, M. Annunziato, and A. Borzì.  Stochastic modelling and control of antibiotic subtilin production.  Journal of mathematical biology, 73(3):727–749, 2016.
  * [32] M. Tosolini, A. Kirilovsky, B. Mlecnik, T. Fredriksen, S. Mauger, G. Bindea, A. Berger, P. Bruneval, W.-H. Fridman, F. Pagès, et al.  Clinical impact of different classes of infiltrating t cytotoxic and helper cells (th1, th2, treg, th17) in patients with colorectal cancer.  Cancer research, 71(4):1263–1271, 2011.