Mathematical modelling of autoimmune myocarditis and the effects of immune checkpoint inhibitors

Autoimmune myocarditis is a rare, but frequently fatal, side effect of immune checkpoint inhibitors (ICIs), a class of cancer therapies. Despite extensive experimental work on the causes, development and progression of this disease, much still remains unknown about the importance of the different immunological pathways involved. We present a mathematical model of autoimmune myocarditis and the effects of ICIs on its development and progression to either resolution or chronic inflammation. From this, we gain a better understanding of the role of immune cells, cytokines and other components of the immune system in driving the cardiotoxicity of ICIs. We parameterise the model using existing data from the literature, and show that qualitative model behaviour is consistent with disease characteristics seen in patients in an ICI-free context. The bifurcation structures of the model show how the presence of ICIs increases the risk of developing autoimmune myocarditis. This predictive modelling approach is a first step towards determining treatment regimens that balance the benefits of treating cancer with the risk of developing autoimmune myocarditis.


Introduction
Immune-related adverse events, and autoimmunity in particular, are an extremely common side-effect of immune checkpoint inhibitors (ICIs), a class of treatments for which approximately 40% of US cancer patients are eligible [1]. Autoimmune side-effects are experienced by 86-96% of cancer patients treated with ICIs, and in 17-59% of patients these side-effects are severe to life-threatening [2,3]. Myocarditis, inflammation of 5 cardiac muscle tissue, is a rare but frequently fatal side effect of ICIs. Of patients treated with ICIs, about 0.1-1% will develop myocarditis, which is a significant increase from the 0.022% incidence rate in the general population [4,5]. Furthermore, myocarditis proves fatal in 25-50% of cases [6,7]. In common with other inflammatory diseases such as rheumatoid arthritis and asthma, the pathology comes from the failure of inflammation to resolve itself, instead becoming a chronic condition [8]. In autoimmune myocarditis, this steady state represents one of the possible outcomes of tumour-immune interactions, i.e. one of the three Es [12].
Similarly, it is reasonable to expect a mathematical model of autoimmune myocarditis to have at least two non-negative steady states, one representing a healthy state, i.e. resolution, (green in Fig. 1) and one 50 representing chronic inflammatory disease, with high numbers of immune cells and dead/damaged cardiomyocytes, or cardiac muscle cells (red in Fig. 1). A third non-negative steady state, a saddle point, separates their basins of attraction (orange in Fig. 1).We propose that the separatrix between the basins of attraction of the healthy and diseased steady states can be seen as a tolerance threshold in the immune system: any challenge to the immune system that stays below this threshold will be tolerated and will not lead to disease, 55 i.e. the system stays in the basin of attraction of the healthy steady state (left panel in Fig. 1). However, if the threshold is crossed, autoimmune myocarditis develops. In a healthy individual, the body is able to tolerate small challenges to the immune system so the separatrix should not lie to close to the healthy steady state. Within this framework, we hypothesize that the impact of ICIs on the tolerance threshold reflects the increased risk of developing myocarditis when a cancer patient is treated with ICIs [4,5]. Thus, the presence 60 of ICIs should move the separatrix closer to the healthy steady state, and an immune challenge that could be controlled before treatment with ICIs, now leads to disease. This is shown in the left panel in Fig. 1 where the trajectories that approached the healthy steady state before (right panel), now lie above the lowered tolerance threshold and thus approach the diseased steady state instead. This hypothesis for the behavioural characteristics of a model of autoimmune myocarditis is tested in this paper.

Aims and outline
We propose a mathematical model of autoimmune myocarditis based on the known biological pathways underlying the disease, and explore whether such a model can produce the behaviours postulated in Fig. 1.
The cell types included in this model represent a minimal set required to describe the mechanisms underlying the development and progression of autoimmune myocarditis. In addition to increasing our understanding of 70 the immunology underlying autoimmune myocarditis, this model takes a first step towards designing patientspecific treatment regimes that balance the benefit of treating cancer with the risk of developing myocarditis.
Section 2 describes the background immunology needed to develop the model. Section 3 presents the model, the steady states and the parameterisation. In Section 4 we show that the presence of ICIs in the system increases the likelihood of a patient developing autoimmune myocarditis. Lastly, in Section 5 we 75 explore how the administration of ICIs might affect the risk of autoimmune myocarditis for different patients by investigating model behaviour in different parameter regimes. We conclude in Section 6 with a brief discussion and proposal for avenues for future research.

Immunology
Here we provide a brief summary of the immunology of autoimmune myocarditis and the mechanisms of 80 action of nivolumab and ipilimumab, two ICIs commonly used in the treatment of cancer. For more detail, the interested reader is referred to the recent review paper [24].
The damage to cardiomyocytes, cardiac muscle cells, that leads to the development of myocarditis can occur spontaneously or it can be triggered by viral infection, cross-reacting T cells, molecular mimicry, or drugs like ICIs, either acting alone or in combination [25]. When cardiomyocytes are damaged or dying 85 their intracellular content can be released. One of the intracellular proteins found in cardiomyocytes is myosin heavy-chain alpha (MyHCα), which is unique to the myocardium and is believed to be one of the major antigens targeted in autoimmune myocarditis [6]. It is a cryptic antigen, meaning the body does not recognize it as "self" and produces an immune response against it [6].
Independent of what triggered the inflammation, sustained, chronic inflammation driven by the release 90 of self-antigens leads to tissue remodeling and fibrosis, which in turn leads to cardiac dysfunction. This is known as dilated cardiomyopathy, the most severe consequence of myocarditis. Fibrosis is an irreversible process, and ultimately, the cardiac dysfunction caused by it will result in heart failure [9]. It is thus key that inflammation is resolved as quickly as possible in order to avoid this fatal outcome.
The immune response can be roughly divided into a fast, non-specific innate immune response and a 95 slowly developing, antigen-specific adaptive immune response [26]. We briefly describe both below.

Innate immune response
The innate immune response relies on tissue-resident sentinel cells, such as macrophages and fibroblasts, and their signaling to recruit other immune cells, such as polymorphonuclear leukocytes (PMNs) [27].
Macrophages are present in high numbers at the inflammatory site in autoimmune myocarditis [28]. Their function depends on their phenotype, which lies on a spectrum between a pro-inflammatory extreme, M1, and an anti-inflammatory extreme, M2. Macrophages maintain high phenotypic plasticity and can change their behaviour in response to environmental cues. M1 macrophages specialise in antigen presentation, phagocytosis of pathogens and dead cells, and stimulate the T cell-mediated immune response, which is further discussed below. M2 macrophages secrete anti-inflammatory cytokines, stimulate wound healing and immune 105 tolerance, and modulate the T helper (Th)2 cell pathway [9,26]. The balance between these two phenotypes needs to be tightly regulated to ensure an effective immune response [29].
PMNs are a group of immune cells that include eosinophils, neutrophils, basophils and mast cells [26].
Eosinophils and neutrophils in particular have been implicated in the pathology of autoimmune myocarditis [30,31]. Both of these immune cell types can respond to antigens by degranulating, a process in which 110 they release toxic compounds into the extracellular space [30,32]. Neutrophils further attack pathogens via phagocytosis or by trapping them in neutrophil extracellular traps [33,34]. Eosinophils play an additional role in regulating the immune response through cytokine secretion [35].
Lastly, fibroblasts are tissue-resident stromal cells that play an important role in maintaining the structural integrity of the tissue. They are essential for wound healing and the development of scar tissue, which 115 they achieve through remodeling of the extracellular matrix [26]. Heart tissue has little-to-no regenerative capabilities so any damaged tissue is replaced by scars, impeding heart function [26,36]. Activation of fibroblasts early on in inflammation can limit the severity of autoimmune myocarditis. If, however, fibroblasts are active for too long, they can cause pathological scarring and tissue fibrosis, which impedes heart muscle function and eventually leads to heart failure [10].

120
The innate immune response is essential in the first few hours to days after tissue injury or infection [26]. It is characterized by high levels of immune-related proteins, PMNs, and edema fluid at the site of inflammation. If the injury or infection cannot be cleared by the innate immune response, the adaptive immune response, driven by antigen-specific T cells, takes over. We summarise this pathway in the next section. 125

Adaptive immune response
Although the adaptive immune response relies on both B and T cells, we focus here on CD4+ T cells as these have been identified as the drivers of the development and progression of autoimmune myocarditis [9].

Antigen presenting cells
Cardiac antigens such as MyHCα that have been released due to cardiomyocyte damage are picked up 130 by antigen presenting cells (APCs). Cells that can function as APCs include macrophages, monocytes and dendritic cells (DCs) [37]. They can be tissue-resident, as is the case for macrophages, or circulate through the body, as is the case for DCs [38]. APCs ingest pathogens, cell debris and other antigens, break them down and present fragments of them on major histocompatibility complexes on their cell surface [26,39].
APCs drain to secondary lymph tissue or the spleen, where they come in contact with mature, naïve T cells [38,40]. These T cells are ready to be activated by an APC that presents to them both the antigen that they recognize, as well as critical co-stimulatory proteins such as CD40, CD80, and CD86 [40]. If a T cell binds an antigen on an APC without simultaneously binding the required co-stimulatory proteins, it goes into anergy, a hyporesponsive state in which the T cell is unable to proliferate or produce cytokines [39,41]. Because the CD4+ T cell-mediated pathway is central to ICI-induced autoimmune myocarditis, we 140 focus here on CD4+ T cells and use "CD4+ T cells" and "T cells" interchangeably.

CD4+ T cells
After activation by an APC, CD4+ T cells go through a process of clonal expansion, during which they divide approximately seven times [42]. The daughter cells produced as part of this process differentiate into a number of phenotypes, most importantly Th and T regulatory (Treg) cells. The phenotype that the T cells 145 assume depends largely on soluble mediators produced by APCs and innate immune cells [39,41,42,43].
There are a number of Th cell subtypes, all of which have a predominantly pro-inflammatory function.
Th1, Th2 and Th17 cells play the most important role in autoimmune myocarditis. Together, these subtypes stimulate the recruitment and activation of macrophages and PMNs, most importantly eosinophils and neutrophils, influence the differentiation of newly activated CD4+ T cells, inhibit each others activity levels, 150 stimulate DC maturation, and damage inflammed tissue, amongst other functions [26,39,44,45].
In contrast, Treg cells are anti-inflammatory T cells that inhibit the proliferation, survival and cytokine secretion of DCs, T cells and innate immune cells, such as PMNs [46,47]. They use a number of different pathways, including cell-cell contact, cytokines and resource competition to control the immune response by influencing the activity of Th cells, DCs, and innate immune cells [48,49]. APCs and T cells thus form an 155 immune pathway targeted at a specific antigen.
Taken together, the innate and adaptive immune responses, and the many links between them, form the basics of the immunology underpinning the development and progression of autoimmune myocarditis. For clarity, much detail has been omitted in this summary, but even so the immense complexity of the underlying mechanisms of autoimmune myocarditis is clear. Fig. 2 provides a graphical overview of the immunology 160 discussed above.

Immune checkpoint inhibitors
The motivation behind ICIs as a cancer therapy is to leverage the powerful immune response against malignancies like tumours. To do this, inhibitory pathways affecting the activation and inhibition of immune cells like CD4+ T cells, are blocked using therapeutic antibodies. Immune-related adverse events arise as 165 side-effects of ICIs when the disinhibited immune response targets non-malignant tissues. Cardiotoxicity, when the immune response is misguided to the heart, is an uncommon irAE, but one with an often sudden onset and quick escalation [7]. The time to onset of symptoms of ICI-induced autoimmune myocarditis varies widely between patients, with some showing symptoms as early as three days after infusion and others taking over a year to display symptoms. The median time to onset is 30-65 days, which corresponds to between one and three doses of ICIs on current treatment schedules [50].
Ipilimumab is one ICI with autoimmune myocarditis as a potential side effect [51,52]. It is an antibody against cytotoxic T-lymphocyte-associated protein 4 (CTLA4). CTLA4 is found on the surface of activated T cells. On Treg cells in particular it is constitutively expressed [46,53]. CTLA4 binds CD80 and CD86, the co-stimulatory proteins on the surface of DCs, with higher affinity than CD28, which is present on the  Another ICI with autoimmune myocarditis as a known side effect is nivolumab, which has a different mechanism of action [51]. Nivolumab is an antibody against a protein called programmed cell death protein-185 1 (PD-1). PD-1 is carried by Th cells on their cell surface [57]. One of its ligands, PD-L1, is found on the surface of peripheral cells, including cardiomyocytes, where it is expressed in response to pro-inflammatory cytokines. A second ligand, PD-L2, is expressed on DCs in response to antigen uptake. Its binding by PD-1 protects the DC against the cytotoxicity of activated T cells and prevents their overactivation [57,58]. When PD-1 binds one of its ligands, the T cell is inhibited and anergy is induced [53]. PD-1-blocking therapy leads 190 to an increase in the expression of proteins for proliferation, activation and effector functions in CD4+ T cells, and thus increases both anti-tumour T cell activity and the probability of inflammation in other organs [59].
Finally, we note that autoimmune myocarditis is often caused by a combination of triggers. Although ICIs increase the risk of developing this disease by increasing T cell activation and decreasing inhibition, 195 they are generally not sufficient to trigger it alone. Often cardiac antigen-specific T cells are in circulation, either due to low negative selection in the thymus or because tumour and gut bacteria-specific T cells are cross-reactive [53,60]. These cells act as secondary trigger for the development of autoimmune myocarditis.
To better direct future research to avenues that are most promising, we need to determine which immunological processes are essential for the development and progression of autoimmune myocarditis, which are of 200 secondary importance, and what roles ICIs play in driving the disease. To aid in making these determinations, we have developed a mathematical model of this disease, which we present in the next section. This model strikes a balance between simplicity, to allow analytical methods to be applied in its study, incorporating cell types determined to be essential by experimental research, and biological realism, to preserve disease characteristics observed in patients.

The model
Here, we introduce the model for autoimmune myocarditis and the effects of ICIs. Our guiding aim is to capture the salient immunological mechanisms and preserve the characteristics of the disease observed in patients while keeping the model tractable. The relative simplicity of ODE models combined with their ability to capture complex, highly nonlinear immunological mechanisms, makes them a suitable starting point 210 for modelling this disease.

The mathematical model
The dynamics of dead/damaged cardiomyocyte numbers are governed by where a 1 , a 2 , R 1 , R 2 and d 1 are positive parameters, and C 0 is the nonnegative initial number of damaged cardiomyocytes. The first term on the right-hand side (RHS) of Eq. (1) where a 4 , a 5 and d 3 are positive parameters, and I 0 is the nonnegative initial number of innate immune cells.

The first term on the RHS of Eq. (2) describes the stimulation of innate immune cells by pathogenic CD4+
T cells [45,64]. Again, pathogenic CD4+ T cells are inhibited by Treg cells. The second term on the RHS 240 of Eq. (2) describes the recruitment of innate immune cells by dead/damaged cardiomyocytes through the innate immune pathway [26,65,66]. The last term describes the death or deactivation of innate immune cells [67,68].
In this model, the pathogenic CD4+ T cell population is assumed to consist primarily of Th cells. The dynamics of pathogenic CD4+ T cells is governed by the ODE where φ and R 3 are positive parameters, P 0 is the nonnegative initial number of pathogenic T cells, a 6 (D N , D I ) is a monotonically increasing function of both nivolumab and ipilimumab doses, and d 4 (D N ) is a monoton-245 ically decreasing function of nivolumab dose. In an ICI-free environment, a 6 (0, 0) > 0 is the activation rate of T cells due to cardiac damage and d 4 (0) > 0 is the deactivation/death rate of pathogenic CD4+ T cells.
The first term on the RHS describes the activation of inactive T cells by dead/damaged cardiomyocytes and their subsequent clonal expansion [26,69]. A fraction 0 < φ < 1 of activated T cells then adopt a pathogenic phenotype [6,26,70]. The processes of activation and differentiation are inhibited by Treg cells, 250 hence the denominator R 3 +R where R 3 is the number of regulatory immune cells required for 50% inhibition of pathogenic CD4+ T cell activation [26,46,48,61]. The second term on the RHS of Eq.
These mechanisms of action are reflected in the model by the effects of D N and D I on the parameters a 6 255 and d 4 .
The number of Treg cells evolves according to the ODE dR dt = a 6 (D N ,

Steady states and linear stability 265
The model described in Section 3.1 was built based on the immunology underlying the progression of autoimmune myocarditis. To assess whether this model can express the behaviour that we hypothesized to be characteristic of a model for this disease (see Section 1. and Fig. 1), we now seek an analytical expression for the steady states of the system. As discussed in Fig. 1 The origin is a steady state of the system (1)- (4). To determine the non-zero steady states of the model, we set the RHSs of Eq. (1)-(4) equal to zero, and use algebraic manipulation to eliminate C, I and R, and establish a polynomial equation satisfied by the non-zero steady states in terms of P . The non-zero steady state numbers of pathogenic CD4+ T cells, P * , satisfy the cubic equation where and The steady state cell counts of the remaining variables are then given by The linear stability of the non-negative steady states is determined by the eigenvalues of the Jacobian matrix

Parameters
Our overarching aim is to develop a model capable of recapitulating the increased risk of developing 275 autoimmune myocarditis upon application of ICIs. The mathematical model has 13 parameters (see Table   1). The grouping of different types of immune cells into classes and the consolidation of many steps in the immunological process into single terms in the model means the model is simple enough to admit some analytical approaches, but it makes it difficult to estimate most parameter values directly from the literature.
Instead, we seek parameter values that give rise to the hypothesized number and types steady states required 280 for the model behaviour illustrated in Fig. 1, in an ICI-free context (D N = D I = 0). The values of parameters d 4 and d 5 can be obtained from literature [73]. For the remaining 11 parameters, we find suitable ranges based on a combination of literature insights and constraints that give rise to the behaviour depicted in Fig.   1. Recall, we proposed that a model of autoimmune myocarditis should have two non-negative, stable steady states and one saddle point, whose position in phase space is influenced by the presence of ICIs. In this 285 section, we explore whether, and for what parameter values, the model can display these characteristics.
We seek parameter sets that satisfy the following criteria (summarised in Table 2): 1. The parameter set must produce three nonnegative steady states. Two must be stable to small perturbations, and one must be a saddle point so that the behaviour of the model captures the ideas in Fig.   1. The origin corresponds to the "healthy" state, with no inflammatory cells. This steady state must be stable. At the origin where [C * , I * , P * , R * ] = [0, 0, 0, 0], the Jacobian matrix is from which we derive the following condition for stability: 2. The saddle point cell counts of damaged cardiomyocytes, innate immune cells and pathogenic CD4+ T cells should be between 10 2 and 10 4 cells, to ensure that a moderately sized perturbation in CD4+ T cells is needed to launch an immune response. Although it is unknown how many immune cells or 290 dead/damaged cardiomyocytes are needed to trigger an immune response, it is known that approximately 3000 DCs are required for a 50% chance of launching an immune response, and we assume that the threshold is at a similar cell count for the other cell types [42,74].
3. The maximum number of damaged cardiomyocytes at any time point must be below 10 6 . The heart consists of approximately 2 − 3 × 10 9 cardiomyocytes in total. To be biologically realistic, the model 295 should not predict the entire heart to be destroyed within a matter of days [75]. [ 42,74] Number of damaged/dead cardiomyocytes should never exceed 10 6 .
[75] Numbers of T cells peak around seven days.
[13] All terms must have significant impact on model output.
[73] To find parameter sets that satisfy criteria one and two, we solve Eq. is employed as a means to perturb the system away from a healthy state as this number of pathogenic T cells will ensure the development of autoimmune myocarditis, as required for testing criteria three and four [53,60]. From these simulations, we extract the maximum number of damaged/dead cardiomyocytes in the system at any point and discard those parameter sets for which this cell count is too high. Next, we select 320 parameter sets for which the number of T cells in the system peaks as close to seven days as possible. The remaining parameter sets thus satisfy criteria three and four. We then select one set of parameter values from the many that satisfy the criteria listed above to make progress on modelling autoimmune myocarditis.
The parameter values in this "default parameter" set are listed in Table 1.
To confirm that the default parameter set satisfies criterion five, we apply a parameter sensitivity analysis  Table 1. In each plot, du denotes the dummy parameter. Results obtained for initial conditions [C 0 , I 0 , P 0 , R 0 ] = [0, 0, 10 4 , 0] cells, at t = 60 days and for 6,000 runs. Asterisk indicates a PRCC value that is significantly different from the dummy parameter, p ≤ 0.01, and the parameter thus has a significant impact on the output variable.
value, and whether or not this impact is significant compared to a dummy parameter. We choose this method because the computational costs are relatively low compared to other methods, such as the extended Fourier  Fig. 4 show that all variables are sensitive to the values of all parameters, and that this parameter set thus satisfies criterion five.

Model simulation with default parameter set
We expect model behaviour to change significantly depending on the magnitude of the perturbation away from the healthy steady state. For small perturbations, the system should be able to recover to the healthy 345 steady state, but for perturbations above a certain threshold the system should tend to the diseased steady state. To explore this, we simulate the model with the default parameter set, setting all initial conditions to zero, except for P 0 which is varied to simulate different immune perturbations away from the healthy steady state at the origin. The time traces in Fig. 5 show that between 250 and 500 pathogenic CD4+ T cells are required to trigger an immune response, which falls within the range of 100-10,000 T cells estimated based on 350 existing literature [42,74]. When autoimmune myocarditis develops, the T cell response takes approximately seven days to launch, while the innate immune response develops almost instantaneously. These behavioural characteristics confirm that the model, simulated using the default parameter set, is capable of reproducing experimental observations [13,26]. Thus, with this particular parameter set, the minimal set of immune cell types and pathways included in this model suffices to recapitulate the disease characteristics seen in patients, 355 and to produce the behaviour required from a mathematical model of autoimmune myocarditis, as postulated in Fig. 1, in an ICI-free context.

The effects of ICIs
In the previous section we established a mathematical model that is rooted in the experimental literature, and can replicate observations of the development of autoimmune myocarditis. We now turn to consider the

Nivolumab
The left panel in Fig. 6 shows what happens to the number of steady states when nivolumab is applied to the system, changing the values of a 6 (D N , D I ) and d 4 (D N ). As previously discussed, nivolumab increases the rate of activation of T cells and decreases the rate of deactivation/death of pathogenic CD4+ T cells.
Thus, as the dose of nivolumab is increased, the value of a 6 (D N , D I ) increases and the value of d 4 (D N ) 375 decreases, moving the system from the default parameter set (starred in the left panel in Fig. 6) towards the bottom right-hand corner. If enough nivolumab is applied, the system crosses a bifurcation which engenders a switch from three steady states to only two steady states. To explore the type of bifurcation that occurs as nivolumab is applied, the values of a 6 (D N , D I ) and origin unstable. Past this bifurcation point, the system only has two steady states: a stable steady state with high immune cell counts and an unstable steady state at the origin. This means that any perturbation to the system away from the healthy steady state leads to the development of autoimmune myocarditis. Applying nivolumab to the system thus initially increases the likelihood of developing autoimmune myocarditis by lowering the threshold for immune perturbations to push the system across the separatrix until, for a high enough dose, any perturbation will lead to disease because the origin becomes unstable.

Ipilimumab
As the dose of ipilimumab is increased, both the activation rate of T cells, a 6 (D N , D I ), and the death rate of Treg cells, d 5 (D I ), increase. The left panel in Fig. 7 demonstrates that the system moves away from the default parameter set (starred) towards the top-right corner as more ipilimumab is applied. As  To explore the type of bifurcation that occurs as the dose of ipilimumab is increased, a 6 (D N , D I ) and d 5 (D I ) are increased simultaneously while the number of steady states and their linear stability are tracked (see right panel in Fig. 7). Again, it is assumed here that the effect of ipilimumab on these parameters is identical in terms of the percentage change in parameter value as a function of ipilimumab dose. The grey dashed line indicates where a 6 (D N , D I ) and d 5 (D I ) attain their default values. As the dose of ipilimumab is increased, the system moves to the right away from the dashed line. Similar to nivolumab, increasing doses of ipilimumab initially lower the threshold for autoimmune myocarditis to develop before the system crosses a transcritical bifurcation and the origin becomes unstable, making the development of the disease almost inevitable. To reach this bifurcation, enough ipilimumab must be added to increase a 6 (D N , D I ) to almost three times its default value.

415
In summary, applying ICIs in increasing doses can thus drastically change the linear stability of the steady states of the model. The presence of increasing levels of nivolumab or ipilimumab leads to a lower threshold for the development of autoimmune myocarditis as the saddle point moves towards the healthy steady state.
For high enough levels of ICIs, a bifurcation occurs, leading to an unstable origin and a system in which the smallest perturbation away from the healthy steady state will lead to disease.

Patient-specific responses
Autoimmune myocarditis is a rare side-effect of treatment with ICIs, so the transition of the system towards resolution or chronic inflammation will differ between patients and per ICI dose. Ultimately, the aim of this research is to be able to predict how the levels of ICIs a patient can receive without developing autoimmune myocarditis, or how different patients might respond to the same treatment protocol. With the 425 parameterised model presented in this paper, we can take preliminary steps towards this goal.   cells are used to simulate a small immune perturbation which, in an ICI-free environment, pushes the system close to the separatrix but not across it. The top-left panel in Fig. 9 shows that in an ICI-free environment the outcome is always resolution, regardless of small changes to the value of a 1 . ICI doses are simulated by a 10% increase or decrease in the appropriate parameter values. This ensures that although the threshold 455 for developing autoimmune myocarditis is lower than in an ICI-free environment, the system has not crossed the transcritical bifurcation that occurs for high doses of ICIs, so the perturbation away from the healthy steady state introduced by the initial conditions is not guaranteed to lead to the development of disease. For the simulation of combination therapy a 6 (D N , D I ) is increased by 21% to take into account the effects of both ICIs, which are assumed to independently affect a 6 . The bottom left panel in Fig. 9 shows that the showing that individual patient characteristics are important to consider when it comes to designing treatment regimen.

Discussion
In this paper, we have presented the first model of autoimmune myocarditis and the effects of ICIs on its development and progression. The model describes how the numbers of damaged/dead cardiomyocytes, We have shown that the addition of ICIs makes it more likely that autoimmune myocarditis develops, i.e. a smaller perturbation away from the healthy state is required to trigger the disease. Furthermore, the number of immune cells at the disease state is higher when ICIs are present, suggesting more severe disease.
We have also shown that small changes in parameter values, representing slight differences between individual 480 patients, can influence whether a particular individual experiences resolution or chronic inflammation. For a certain ICI dose or immunological perturbation, the model predicts that some patients may be pushed towards chronic inflammation while others return to a healthy state.
Our model predicts that for a large enough dose of either nivolumab or ipilimumab, the development of autoimmune myocarditis is inevitable for even the smallest perturbation away from the healthy steady state. of the perturbation, i.e. the number of active pathogenic CD4+ T cells added to the system, can be carefully controlled. Results of such an experiment could be used to validate the mechanisms proposed in this paper.
We highlight that the parameter values relevant for this model are highly uncertain. No data with which to calibrate the model is currently available, and the grouping of many cell types and interactions into single 490 variables and terms makes it difficult to estimate parameter values from existing literature. As the magnitude of the perturbation away from the healthy steady state significantly impacts the timescale on which the cell counts evolve, and the timescale on which pathogenic T cells are activated is one of the criteria around which the model is parameterized, the parameter values resulting from the method we have employed are clearly dependent on the size of the perturbation. Any change to the assumed magnitude of the perturbation may 495 require changes to the parameter values to ensure that all criteria are still met. However, we highlight that the qualitative behaviours of the system are not affected.
Should patient-specific data with which to calibrate become available, this model could be used to estimate the maximum tolerated dose of ICIs. Fitting the model parameters to patient-specific data allows us to compute where the saddle point lies, how far away it is from the healthy steady state and thus how big a 500 perturbation away from the healthy steady state a patient can tolerate. This enables assessment of the risk of a patient developing autoimmune myocarditis for a certain dose of ICIs. It also facilitates predictions about outcomes for individual patients and the timeline on which disease may develop. This is all useful data that could potentially inform the treatment plan for specific patients.
To be able to have confidence in model predictions, we must ensure that the parameters of the model are 505 identifiable given the available data. To obtain an understanding of the type and quality of data required to establish parameter values for this model, we conducted a parameter identifiability analysis. Structural identifiability analysis performed using DAISY indicated that the model is only structurally identifiable if all four variables are measured as outputs. To assess the practical identifiability of the model parameters, the maximum loglikelihood profiles of subsets of the 13 parameters in the ICI-free model are computed, using 510 the methods explained by Simpson et al. [77]. An example of the results of the practical identifiability analysis are shown in Fig. 10, where the profile loglikelihoods show that the death rates and φ form a practically identifiable subset of parameter space. More details regarding the methods and results of the practical identifiability analysis can be found in Appendix A.
In the future, this model may be used together with methods from optimal control theory to predict 515 optimal drug regimen. If the goal is to maximally attack the cancer with ICIs while preventing the development of autoimmune myocarditis, the frequency and quantity of ICI dosage needs to be carefully planned.
Our results demonstrate that as the drug dose increases, the saddle point approaches the healthy steady state, which leaves less and less room for fluctuations in the number of DCs and pathogenic immune cells without crossing over into the basin of attraction of the disease steady state. Optimal control theory can be 520 used to determine where the ideal balance between the opposing aims of treating the tumour and preventing autoimmunity lies. With this model, we have taken the first steps in consolidating the current knowledge of autoimmune myocarditis and the effects of ICIs on its development and progression. It will be important to further investigate which cell types and molecules are the critical players in this disease in order to work towards 525 prevention of this side effect of ICIs. Though much work remains to be done, we have started on the path towards a mathematical model of autoimmune myocarditis and the role of ICIs that can be used to determine which cell types, cytokines and other components are critical in the development and progression of this severe side effect. With these essential immune components identified, work can start towards reducing the risk of severe cardiotoxic side effects of ICIs in patients.

Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The large number of parameters and relatively small number of variables in the model makes a comprehensive practical identifiability analysis infeasible. The parameter set is therefore split into subsets to investigate which subsets of model parameters are identifiable. We use the profile loglikelihood method to explore the practical identifiability of parameters [77]. In this method, the parameter of interested is varied over a set range and for each value, the values of the other parameters in the subset are optimized to maximize the 545 loglikelihood.

S.A. van der
A synthetic data set is obtained by simulating the model for 60 days, obtaining an "observation" of all four variables on days 0, 1, 2, 5, 10, 30 and 60. Structural identifiability analysis performed with DAISY indicates that the model parameters are only structurally identifiable when all variables are observed, therefore observing only a subset of the variables will guarantee practical unidentifiability. The observations are 550 concentrated in the first few days so the fast early dynamics of the system are captured. Normally distributed measurement noise with a standard deviation of 5% of the cell count observed is added to all data points. where y j is the model output of the variable j for a given parameter set, y d is the data set, i denotes one of n data points, and SD is the standard deviation for a specific data point, computed as SD = 0.05y d i,j . The results in Fig. 10 show that the death rates and φ form a practically identifiable subset with the 560 synthetic data set created with the above requirements. The set of the activation rates, however, do not form a practically identifiable subset given the specific settings used here, as is shown in Fig. A.12 where a subset of the activation rates consisting of the parameters a 1 , a 2 , a 4 , and a 5 is shown to be practically unidentifiable.
The subsets can also be formed by grouping parameters by ODE they appear in. As is to be expected, this  leads to practically unidentifiable subsets as parameter values in these subsets can be adjusted to compensate for the varying value of the parameter of interest. An example of this is shown in Fig. A.13, which shows the flat profile loglikelihoods for the parameter subset a 1 , a 2 and d 1 , all of which appear in the ODE describing the dynamics of damaged/dead cardiomyocytes.