key: cord-0020254-635amtra
authors: Le, Trang; Su, Sumeyye; Shahriyari, Leili
title: Investigating Optimal Chemotherapy Options for Osteosarcoma Patients through a Mathematical Model
date: 2021-08-06
journal: Cells
DOI: 10.3390/cells10082009
sha: 6b76f573c1bead74f26f49a6b3a7f6676f639fff
doc_id: 20254
cord_uid: 635amtra

SIMPLE SUMMARY: Osteosarcoma is a rare type of cancer with poor prognoses. However, to the best of our knowledge, there are no mathematical models that study the impact of chemotherapy treatments on the osteosarcoma microenvironment. In this study, we developed a data driven mathematical model to analyze the dynamics of the important players in three groups of osteosarcoma tumors with distinct immune patterns in the presence of the most common chemotherapy drugs. The results indicate that the treatments’ start times and optimal dosages depend on the unique growth rate of the tumor, which implies the necessity of personalized medicine. Furthermore, the developed model can be extended by others to build models that can recommend individual-specific optimal dosages. ABSTRACT: Since all tumors are unique, they may respond differently to the same treatments. Therefore, it is necessary to study their characteristics individually to find their best treatment options. We built a mathematical model for the interactions between the most common chemotherapy drugs and the osteosarcoma microenvironments of three clusters of tumors with unique immune profiles. We then investigated the effects of chemotherapy with different treatment regimens and various treatment start times on the behaviors of immune and cancer cells in each cluster. Saliently, we suggest the optimal drug dosages for the tumors in each cluster. The results show that abundances of dendritic cells and HMGB1 increase when drugs are given and decrease when drugs are absent. Populations of helper T cells, cytotoxic cells, and IFN- [Formula: see text] grow, and populations of cancer cells and other immune cells shrink during treatment. According to the model, the MAP regimen does a good job at killing cancer, and is more effective than doxorubicin and cisplatin combined or methotrexate alone. The results also indicate that it is important to consider the tumor’s unique growth rate when deciding the treatment details, as fast growing tumors need early treatment start times and high dosages.

Osteosarcoma is a rare kind of cancer that occurs in the bone, and around 1000 new cases of osteosarcoma are identified each year in the United States [1] . It can affect people of any age, but it mostly occurs in children who are 10 to 14 and in adults who are 65 and older [2, 3] . There are some factors such as gender, age, heritable syndromes, and certain other conditions that affect the risk of osteosarcoma [4] . The types of standard treatment for osteosarcoma include surgery, chemotherapy, radiotherapy, and targeted therapy [5] .

Although neoadjuvant chemotherapy has shown improved results in treating osteosarcoma, patients with metastases have continued to have low survival rates [6] [7] [8] . Immunotherapy and targeted therapy are known alternative treatments of osteosarcoma; however, they are still inefficient for many patients [9] . Meanwhile, resistance to radiotherapy has also been observed in osteosarcoma tumors [10, 11] . To overcome therapeutic growth, that is, the rate of change of cancer cell population is proportional to [C] 1 − [C] C 0 , with C 0 being the carrying capacity. Cytokines in [58] were modeled to be proportional to the populations of cells that produce them.

The full system of ODEs from the model in [58] is:

Here, λ parameters denote proliferation, activation, or production rates. δ parameters denote inhibition, decay, or death rates. A M N , A T N , and A D N correspond to the production/proliferation rates of naive macrophages, naive T cells, and naive dendritic cells, respectively; and α NC corresponds to the fraction of dying cancer cells that become necrotic cells. The descriptions of all variables in this system are given in Table 1 .

We build upon this system of ODE by adding the interactions of the variables in this system with the following chemotherapy drugs: methotrexate, doxorubicin, and cisplatin. The interaction network of these drugs with cells and cytokines of osteosarcoma tumor microenvironment is shown in Figure 1 . We used an exponential kill model, as introduced in [84] , to describe how chemotherapy affects the cancer microenvironment, and modeled the change in population of the new model's variables over time in the unit of days. The details of the effects of chemotherapy drugs on immune cells and cancer cells are explained below (changes to Equations (1)- (14) are in bold). 

Name Description

Helper T-cells T C

Regulatory T-cells D n Naive dendritic cells D Activated dendritic cells antigen presenting cells M N Naive macrophages includes naive macrophages and monocytes M Macrophages includes M1 macrophages and M2 macrophages C Cancer cells N Nectrotic cells H HMGB1 µ 1 Cytokines group µ 1 includes effects of TGF-β, IL-4, IL-10 and IL-13 µ 2 Cytokines group µ 2 includes effects of IL-6 and IL-17 I γ IFN-γ A 1 methotrexate methotrexate concentration at tumor site 

All chemotherapy drugs in our model aim to kill tumor cells, though they have different mechanisms of action. Methotrexate hinders DNA synthesis in fast dividing cancer cells by inhibiting folate-dependent pathways [85] . Doxorubicin can kill cancer cells by binding to DNA-associated enzymes, intercalating the base pair of DNA's double helix, and targeting many molecular targets such as topoisomerase enzymes I and II, which results in DNA damage [86] . Cisplatin binds platinum to DNA by forming inter-stranded and intra-stranded crosslinks, and thus induces DNA damage which leads to cell death in rapidly proliferating cells [87, 88] .

Similarly to [84, 89] , we use saturation kill term (1 − e βA ) to model the direct cytotoxic effects of chemotherapy drugs on cancer cells, where β is the drug efficacy parameter, and A is the drug concentration at the tumor site. This is based on the observation that at low concentrations, the cancer killing effects of these drugs are almost linear, but at very high concentrations the cancer killing effects plateau. Unlike doxorubicin and cisplatin, methotrexate can only eliminate cancer cells during certain phases of the cell cycle, so we added the term ( f − τ a + 1 24a ) to methotrexate's cytotoxic effect to account for this phenomenon, as modeled in [84] . Here, f denotes the fraction of cells in the vulnerable phase of the cell cycle for methotrexate, a denotes cell cycle time in days, and τ is defined to be minimum(T, f a), with T being drug exposure time in days.

Besides its direct role in targeting tumor cells, cisplatin has also been reported to increase cytotoxic cells' cancer killing capability by upregulating MHC-1 expression in cancer cells [35, 87, 90, 91] . We also use a saturation term to describe this effect, as it is very likely that a high concentration of cisplatin can also plateaus in the upregulation of MHC-1 in cancer cells. We make the assumption that the concentration of cisplatin at which this effect slows down is about the same concentration at which the cancer killing effect of cisplatin slows down, so we use the same drug efficacy parameter β 3 in both terms, resulting in the following equation for cancer cells:

where δ CT c A 3 represents cisplatin's promotion of cytotoxic cells' cancer killing ability; K C is rate of chemotherapy-induced tumor death; and β 1 , β 2 , and β 3 are the medicine efficacy coefficients of methotrexate, doxorubicin, and cisplatin, respectively. A description of every chemotherapy-related parameter in our model is given in Table 2 .

As a proportion of cancer cells killed by chemotherapy drugs become necrotic cells, we describe the change in population of necrotic cells with the presence of chemotherapy as follows:

where α NCA is the fraction of dying cancer cells induced by chemotherapy that turn into necrotic cells.

Since chemotherapy does not only eliminate tumor cells but also kills immune cells, we include the effects of chemotherapy in the equations of immune cells as well. Similarly to [89] , we assume that the same quantity of chemotherapy drugs is required to affect cancer cells and immune cells, even when the rate at which chemotherapy kills cancer cells is different than the rate at which it kills immune cells. Hence, we use the same drug efficacy coefficients for cancer and immune cells, but different rates of drug-induced cell death, leading to the following modified immune cell equations:

(1 − e −β 1 A 1 )

where K M N , K M , K T N , K T h , K T r , K T c , K D N , and K D are the rates of chemotherapy-induced cell death of naive macrophages, macrophages, naive T cells, helper T cells, regulatory T cells, cytotoxic cells, naive dendritic cells, and dendritic cells, respectively.

Chemotherapy drugs are given through IV infusion in osteosarcoma treatments, so their bioavailability is 100%. Thus, we use the following equations to model the changes in concentration of the chemotherapy drugs at the tumor site over time:

Here, v A 1 (t), v A 2 (t), and v A 3 (t) are the amounts of methotrexate, doxorubicin, and cisplatin injected per day per liter of body volume, with the unit of mg/L per day; δ A 1 , δ A 2 , and δ A 3 are the decay rates of methotrexate, doxorubicin, and cisplatin, respectively.

A recent study showed that among the available digital cytometry methods with which to deconvolve the gene expression data of a tissue into relative abundances of cells in the tissue, CIBERSORTx B-mode performs best [92] . In our previous study [58] , we used the gene expression data of 88 osteosarcoma tumors from TARGET data set and applied CIBERSORTx B-mode to estimate the fractions of immune cells in each of these tumors. Then, k-means clustering was performed on the estimated cell fractions, and three clusters of osteosarcoma tumors with distinctive immune patterns were found. In this work, we study the impacts of chemotherapy drugs on the same three clusters. The average immune fractions in each cluster are presented in Figure 2 . To obtain immune cell populations, we multiplied the immune fractions estimated from CIBERSORTx by a scaling factor α dim . From the supplementary data of TARGET project, which include information on the percentage of tumor, necrotic, stroma, and normal cells of each sample, we derived the cancer and necrotic populations from immune populations. We used the percentage of normal cells to denote the percentage of total immune cells. In particular, given total immune population I, we computed the populations of cancer and necrotic cells as follows:

where C and N are the populations of cancer and necrotic cells, respectively.

We chose α dim based on the average cancer cell population of osteosarcoma. The mean volume of Ewing sarcomas is reported to be 275 mL [93] , and a study found Ewing sarcoma volumes to be similar to osteosarcoma volumes [94] . Therefore, we approximated the mean volume of osteosarcoma to be 275 mL. Since osteoblasts have a diameter of 20-50 µm [95] , we estimated that osteosarcoma cells, which are malignant osteoblasts, have an average diameter of 35 µm, leading to an estimated average of 6.4 × 10 9 cancer cells in osteosarcoma tumors. Thus, α dim was chosen to be 1.765 × 10 8 to match the average cancer cell population among all patients in TARGET data to 6.4 × 10 9 . It is worth noting that a smaller size of osteoblasts has also been reported [96] . However, α dim is simply a scaling factor and has no impacts on the dynamics of cells and cytokines in our system; thus, the size of osteoblasts does not affect our results.

Given a treatment regimen of interest, we applied its standard dosage to our model. Most doses of chemotherapy drugs for osteosarcoma are measured in mg/m 2 , but we modeled the drug concentration at the tumor site in milligrams per liter of body volume. We therefore need to convert drug doses from mg/m 2 to mg/L. We used an average body surface area of a human male of 1.9 m 2 [97] and an average male body volume of 59.7 L [98] for conversion. That is, for example, 75 mg/m 2 would be equivalent to:

The drug efficacy coefficients, cell cycle time, and fraction of cells in the vulnerable phase of the cell cycle are taken from [84] . Using the molecular mass of chemotherapy drugs [89, [99] [100] [101] , we can convert the drug efficacy coefficients given in [84] to units of (mg/L) −1 : 

Drug efficacy coefficients for doxorubicin-resistant and cisplatin-resistant cells were also included in [84] , and can be converted in a similar way. The values for all chemotherapyrelated parameters in our model and their sources are given in Table 2 .

The fraction of cancer cells killed by chemotherapy, K C , was taken from [89, 102] to be 0.9, based on the notion that chemotherapy strength is one log kill [103] . Since chemotherapy is more effective at eliminating fast proliferating cells, it is safe to assume that the rates of chemotherapy-induced death of immune cells are smaller than those of chemotherapy-induced death of cancer cells. Hence, K M N , K M , K T N , K T h , K T r , K T c , K D N , and K D are assumed to be smaller but in the same order of magnitude as K C , and we use a value of 0.6 for all of them, as in [102] . Decay rates of chemotherapy drugs were derived from their elimination half lives in the following way:

where δ A is the decay rate of A. On average, the elimination half lives of doxorubicin, cisplatin, and high-dose methotrexate are 2 hours, 25 min, and 11.5 h respectively [104] [105] [106] , resulting in the corresponding decay rates of 8.3178, 39.9253, and 1.4466.

As no values for α NCA and δ CT c A 3 can be found in the literature, we assume biologically reasonable values for these parameters. α NCA is the fraction of dying cancer cells induced by chemotherapy that become necrotic cells, so it is bounded between 0 and 1. We make the assumption that a large proportion of dying cancer cells from treatment turn into necrotic cells, and set α NCA = 0.8. For δ CT c A 3 , we assume that cisplatin at maximum effect can double the cancer killing ability of cytotoxic cells, or equivalently δ CT c A 3 = 1. In order to investigate whether our assumptions for those parameters have large impacts on the cancer cell population, we performed sensitivity analysis and studied the change in cancer cell population after treatment while varying these two parameters. 

To achieve additional numerical stability, non-dimensionalization of the whole system was carried out. For each variable X of the original system in [58] , its dimensionless form can be written as:

where X ∞ is the steady state value of X given in [58] . For the newly added variables, which are the chemotherapy drugs, we introduced the following non-dimensional variables:

where A is the dimensional variable, δ A is the decay rate of A, and v * A is the amount of drug A injected on its first injection day of the treatment. Further details on nondimensionalization are given in Appendix A.

To solve the non-dimensional system of ordinary differential equations, we used the solve_ivp function from Scipy package in python [107] , with initial conditions from a chosen data point of interest in each cluster.

We performed local gradient-based sensitivity analysis on all chemotherapy-related parameters to study their impacts on the outputs of the system. For the non-dimensional system dX dt = F(X, θ, t) with model parameters θ = θ 1 , ..., θ N , the local (first order) sensitivity of parameter θ i with respect to the variable X is defined as [108] :

As we are mainly interested in how drug-related parameters affect the number of cancer cells, we calculated the sensitivity of treatment parameters with respect to cancer and total cell populations. Since the effects of the treatment do not reach steady state, we consider time-dependent sensitivity. That is, we measure sensitivity of parameters in every time step throughout the treatment and some time after. The change in sensitivity of θ i over time can be derived as follows:

By applying the chain rule, we have:

In addition, we also look at the relative sensitivity, which is commonly used in metabolic control analysis of biological networks [108] :

Then, we compute the average sensitivity of each type over a period of time T:

The sensitivity varies for different values of the parameters, so we consider a small neighborhood Ω(θ) of the given parameter set and calculate:

where the integrals are computed using a numerical technique called sparse grid points [109, 110] .

In particular, to evaluate S i , sparse grid point method samples θ k from parameter space Ω for k = 1, ..., m, where the number of samples m are chosen by the user. Then S i (θ k ) is evaluated for k = 1, ..., m, and the summation m ∑ k=1 S i (θ k ) is used as an approximation to S i .

Here, to generate sparse grid points samples, we apply Matlab package nwspgr provided by [109] , and we use m = 39.

We introduce a framework to find the appropriate dose of a given treatment regimen for each patient. With it, one aims to achieve a target cancer cell population at a given time t, typically to reduce the size of tumor for the day of surgery. Hence, we want to minimize the difference between the estimated cancer cell population at t from our model and the target cancer cell population. The least square error is used to model this difference, as least square error penalizes high errors more than absolute error, and is less computationally expensive than hinge loss.

Since high doses of chemotherapy are known to induce high toxicity to the patient, we would like to achieve cancer cell populations close to the target population with the smallest dosages possible. Therefore, we added a L1 norm regularization term, which is the sum of absolute values of drug dosages, to the loss function. We used the L1 norm instead of the L2 norm, because the L2 norm penalizes high doses of individual drugs more severely. In most osteosarcoma treatments, methotrexate is used at a much higher dose than other drugs, so L2 norm would put too much weight on methotrexate dose and neglect other drugs' doses.

As patients can only tolerate certain dosages of chemotherapy and drug doses cannot be negative, we also put lower and upper bounds on the drug doses. Thus, we have the following loss function:

Here, v is a vector of length M, denoting the doses of the M drugs in the given treatment. t is the time of interest at which cancer cell population is evaluated for optimization. C(v, t) is the cancer cell population with drug doses v at time t of interest, and is computed via our ODE solver. C target is the target cancer cell population at time t, and is chosen by the user. U i is the highest allowable value for v i . κ is the regularization parameter; the higher κ is, the more the optimizer focuses on achieving small doses and less on achieving small error betweenĈ(v, t) and C target .

The optimize.minimize function from Scipy package in python is used to solve this optimization problem, with the outputs being the optimal doses.

Typical treatments for osteosarcoma include neoadjuvant chemotherapy, usually for 10 weeks, then surgery, and adjuvant chemotherapy after the surgery for up to a year [40] . The most common chemotherapy regimen for osteosarcoma in children and young adults is the MAP regimen, which is a combination of doxorubicin, cisplatin, and high-dose methotrexate [39] . This regimen consists of six 35-day cycles; two cycles are applied before surgery and the remaining four are applied after surgery. In each cycle, 37.5 mg/m 2 of doxorubicin and 60 mg/m 2 of cisplatin are administered through IV per day on days 1 and 2, and 12,000 mg/m 2 of methotrexate is administered through IV over 4 h per day on days 22 and 29 [111, 112] . Different infusion schedules have been used for doxorubicin and cisplatin: doxorubicin can be injected as a bolus or a 4-h infusion each day, or a continuous infusion over 48 h, whereas cisplatin can be injected over 2 or 4 h each day, or continuously over 72 h [112] . We studied the dynamics of cell and cytokine populations in large osteosarcoma tumors during neoadjuvant MAP treatment, which includes two 35-day cycles. In particular, we used the steady state values of cells and cytokines in [58] as initial conditions, the typical dosage of the MAP regimen as the drug inputs, and 4-h infusions on previously specified days as the injection schedule for each drug. We set the start of chemotherapy treatment to be 7 days after biopsy, as it usually takes a few days to receive the results of the biopsy. Figure 3 shows that for all clusters, cancer cell populations are reduced significantly after two cycles of MAP treatment. It is important to note that the cancer cell populations do not reach zero after chemotherapy, so cancer cells will start growing again after chemotherapy. However, the goal of neoadjuvant therapy is not to eradicate cancer cells completely, but only to reduce the boundaries of the tumor and to remove any small metastases that have not been detected [113] . Cluster 2 had the highest cancer cell population at the start of treatment, so naturally cluster 2 also had the highest cancer cell population left after neoadjuvant therapy. Interestingly, cluster 1 had approximately the same number of cancer cells as cluster 3 at the start of treatment, but ended up with a higher cancer cell population than cluster 3 after treatment. This is because in each chemotherapy cycle, there are a few weeks where no chemotherapy drugs are administered in order to allow the patient to recover from the drugs' toxicity, and thus during these few weeks, cancer cells can start growing again. Cluster 1's cancer cell population, which was reported in [58] to have the highest growth rate in the three clusters, grew more during the weeks with no drugs given, resulting in a higher number of cancer cells after treatment compared to cluster 3. This observation suggests that we should take the patient's cancer growth rate into account when choosing their chemotherapy dosage.

During the MAP treatment, necrotic cells, dendritic cells, and HMGB1 oscillated between increasing and decreasing in abundance. Since chemotherapy drugs aim to kill cancer cells, and a fraction of dying cancer cells become necrotic cells, the population of necrotic cells increases on the days the drugs are injected and there are many druginduced dying cancer cells. However, during the weeks where no chemotherapy drugs are administered, only drug-free dying cancer cells can become necrotic cells, and with the cancer cell population being already reduced by the previously given drugs, the number of drug-free dying cancer cells is small, leading to a decrease in necrotic cell population. HMGB1 is mainly produced by necrotic cells, so HMGB1 abundance increases when the necrotic population increases and decreases when the necrotic population decreases.

Meanwhile, dendritic cells are activated largely by HMGB1, so the dynamics of dendritic cells share the same trend as the dynamics of HMGB1. That means both HMGB1 and dendritic cells increase on the days chemotherapy drugs are administered and decrease on the weeks with no drugs given. An increase in dendritic cell population right after doxorubicin [114] [115] [116] [117] or cisplatin [35, 36, 87] administration, and a rise in HMGB1 production following doxorubicin [116] [117] [118] , have been shown in several studies, which aligns with our results.

We observed that in general, helper T cells, cytotoxic cells, and IFN-γ increase in population during chemotherapy. There are many studies that report an increase in helper T cells' and/or cytotoxic T cells' abundance due to doxorubicin [117] [118] [119] [120] [121] [122] , cisplatin [35, 87, 90, 123] or methotrexate [124] , and thus they support our findings. Doxorubicin has been known to induce immunogenic cell death, which leads to the maturation of dendritic cells and accordingly the activation of helper and cytotoxic T cells [114, 120, 125] . IFN-γ is produced by helper T cells and cytotoxic cells; thus, IFN-γ abundance also increases as these two cells increase in population during MAP treatment. The increase in IFN-γ level after administration of doxorubicin and cisplatin has also been observed in multiple studies [87, 116, 117, 126] .

On the other hand, populations of macrophages, regulatory T cells, cytokines groups µ 1 , and µ 2 mainly decrease during MAP treatment. These immune cells are not affected by the necrotic cell death process caused by chemotherapy, so they decrease in population during chemotherapy, as they are also killed by the drugs. µ 1 and µ 2 are produced by helper T cells, macrophages, and cancer cells. Even though the helper T cell population increases during treatment, macrophage and cancer cell populations decrease more so, which leads to an overall decrease in µ 1 and µ 2 throughout MAP treatment. Several other studies have also reported a reduction in regulatory T cell quantity due to cisplatin [35, 87, 90 ] and a decrease in the level of IL-6, which is the main component of µ 2 , due to methotrexate and doxorubicin [124, [126] [127] [128] .

To study the impacts of the newly introduced parameters on the outputs of the model, we performed local sensitivity analysis on the chemotherapy-related parameters using the non-dimensional system. The initial conditions for sensitivity analysis are the large tumors in each cluster, which we used without-treatment steady state values to represent. It is worth noting that the cell cycle time, a, was not included in this sensitivity analysis, because it is a simple measurement rather than a parameter that needs to be estimated or fitted to the experimental data. The most sensitive time-averaged parameters in terms of sensitivity and relative sensitivity are presented in Figure 4 .

In all clusters, the initial fraction of cells in the vulnerable phase of the cell cycle f has the largest impact on cancer cell population among treatment-related parameters according to both the sensitivity and relative sensitivity analyses. The rate of chemotherapy-induced cancer cell death, K C , and the drug efficacy coefficients of doxorubicin and cisplatin, β 2 and β 3 , are also sensitive to cancer and total cell population during treatment. Meanwhile, the drug efficacy coefficient of methotrexate, β 1 , does not seem as sensitive, but the decay rate of methotrexate is.

We notice that the parameters whose values are assumed, α NCA and δ CT c A 3 , do not have large effects on the cancer cell population or total cell population based on the results of the sensitivity analysis. To further confirm this, we also plotted the cancer cell populations after treatment with different values of these parameters. We chose α NCA ranging from 0.2 to 1, because it is a fraction and thus is bounded between 0 and 1, and δ CT c A 3 ranges from 0.2 to 5 times its original value. Figure 4C ,D shows that varying either of these parameters results in negligible changes to cancer cell population after treatment. 

The effectiveness of chemotherapy is highly dependent on the existence of resistant cancer cells. We are interested in studying the changes in quantity of cells and cytokines when osteosarcoma cells are resistant to one or multiple drugs within the MAP regimen. As mentioned in Section 2, we obtained drug efficacy coefficients from [84] ; these values were estimated to fit the survival data of cancer cells under different chemotherapy drugs. The same study [84] also included the estimated drug efficacy coefficients of doxorubicin and cisplatin in doxorubicin-resistant and cisplatin-resistant cancer cells, respectively. Using these parameter values, we plotted the dynamics in the osteosarcoma microenvironment during MAP treatment when cancer cells are resistant to either doxorubicin or cisplatin, or to both drugs. Since methotrexate-resistant cells were not used in [84] , and hence no parameter values were available for them, we did not model the dynamics with methotrexate-resistant cells. Figure 5A shows that MAP treatment is not as effective at shrinking the tumor when cancer cells are resistant to doxorubicin; the cancer cell population after treatment is about 60% to 70% higher in doxorubicin-resistant cells than in non-doxorubicin-resistant cells ( Table 3 ). The smaller reduction in cancer cell population of doxorubicin-resistant cells during doxorubicin administration means fewer necrotic cells are produced in the process, and accordingly a lower level of dendritic cells ( Figure 5A ), as necrotic cells indirectly promote dendritic cell maturation through the release of HMGB1. We notice no clear difference in the dynamics of T cells and macrophages compared to the microenvironment of non-doxorubicin-resistant cells. It is worth noting that we modeled only cancer cells to be resistant to chemotherapy drugs, so immune cells were by design not resistant to these drugs. On the other hand, with cisplatin-resistant cells, we observed little difference in the reduction of cancer cell population compared to non-cisplatin-resistant cells ( Figure 5B , Table 3 ). This is due to the fact that cisplatin's drug efficacy parameter, β 3 , is small compared to methotrexate and doxorubicin's drug efficacy parameters, resulting in cisplatin having a rather minor effect on cancer reduction in the MAP treatment. Hence, the resistance to doxorubicin matters more than the resistance to cisplatin. Since the cisplatinresistance does not have a large impact on the effectiveness of MAP treatment, cancer cells that are resistant to both doxorubicin and cisplatin have similar dynamics to doxorubicinresistant cancer cells ( Figure 5A,C) .

We studied the effects of delays in the starting time of treatments on the tumors' responses to the treatments. Since the tumor's growth rate depends on the tumor's size, we investigated the effects of delaying the chemotherapy treatments in small, medium, and large tumors, separately. Small tumors were chosen as follows: we first chose the tumor with the smallest cancer cell population in cluster 1, then found the tumors in cluster 2 and 3 that had cancer cell populations closest to the chosen tumor from cluster 1. Medium tumors were taken to be the mean values of all patients in each cluster. For large tumors, we took the without-treatment steady state values. We plotted the dynamics of the cancer cell population in each cluster when chemotherapy was started 1 week after diagnosis, which we assume is the earliest start time, as it takes a few days to obtain biopsy results; and the dynamics of beginning 1 month, 3 months, and 6 months after the initial diagnosis. Figure 6 and Table 4 show that in small and medium tumors, the cancer cell populations after treatment are higher the longer we wait to start the chemotherapy. Thus, starting chemotherapy earlier leads to better outcomes with these tumors. On the contrary, the cancer cell population stays the same after treatment in large tumors regardless of the treatment start time. This is because these large tumors are at steady state and do not grow more while the patient waits for the treatment. Theoretically, the treatment start time does not matter as much for tumors at steady state or close to reaching steady state. However, in reality, when tumors are large, the functionality of the cancerous body part is likely compromised, and the quality of the patient's life is affected, which makes us want to start the chemotherapy promptly for large tumors. It was previously observed in [58] that tumors in cluster 1 grow fast even when the tumor is small, whereas tumors in clusters 2 and 3 start growing fast when the tumor is a bit bigger. Hence, when we delay the treatment for a long time, the small tumor in cluster 1 grows quickly and ends up with a far higher cancer cell population after the treatment than in other clusters, as seen in the treatments starting at 3 months and 6 months ( Figure 6A , Table 4 ). For small tumors in clusters 2 and 3, even though they do grow while the patients wait for treatment, their growth rates are not as high, and their cancer cell populations after the treatment are still relatively small despite the treatment delays. Therefore, it is important to start the treatment early for small tumors of cluster 1, and it would be ideal but not as urgent to start the treatment early for small tumors in clusters 2 and 3. 1.34 × 10 10 1.87 × 10 9 1.87 × 10 9 1.87 × 10 9 1.87 × 10 9 Figure 6B indicates that medium-sized tumors in all three clusters grow comparably fast, and since their cancer cell populations after treatment will not be very small, we should start chemotherapy for them as early as possible. We also notice that for small and medium tumors in all clusters, the differences in cancer cell population after any treatment were not significant when comparing starting the treatment 1 week or 1 month after the diagnosis. However, treatments starting 3 or 6 months after resulted in much larger cancer cell populations after the treatment. Based on our model, it is thus not recommended to start the chemotherapy several months after the diagnosis, but rather to start within a month of the initial diagnosis.

We investigated the effects of two other chemotherapy regimens on the osteosarcoma microenvironment. A combination of doxorubicin and cisplatin (AP) is a very common treatment of osteosarcoma tumors in older adults, as they are less likely to be able to tolerate high-dose methotrexate. This regimen consists of three preoperative 21-day cycles, where 25 mg/m 2 of doxorubicin is given as a bolus once per day from day 1 to day 3, and 100 mg/m 2 of cisplatin is given as a continuous infusion over 24 h on day 1 in each cycle [129, 130] . High-dose methotrexate (MTX) has also been used as a single agent to treat osteosarcoma, with four courses of 8 to 12 mg/m 2 given weekly before surgery [131] . In this study, we used the average dose, which is 10 mg/m 2 of methotrexate injected over 4 h on day 1 every week for this regimen. Figure 7 shows that MTX and AP regimens both had higher cancer cell populations after treatment than MAP. This agrees with the finding from [132] : that the AP regimen is less effective but safer than the MAP regimen. Meanwhile, MTX as a single agent was reported to be insufficient as a neoadjuvant therapy for osteosarcoma in [131] , which used the same MTX dosages and schedules as this study. Overall, according to our model, MAP is the superior treatment to MTX and AP in terms of cancer-killing ability. In fact, a recent study reports that MAP is still the favorable option for osteosarcoma among various combinations of chemotherapy drugs [133] .

The AP regimen has relatively similar dynamics of cells and cytokines as the MAP regimen. That is, the populations of HMGB1, necrotic cells, and dendritic cells increase when drugs are given and decrease when no drugs are given; populations of helper T cells, cytotoxic cells, and IFN-γ decrease less than they increase, so in general they increase during treatment; and regulatory T cells, macrophages, and cytokine groups µ 1 and µ 2 generally decrease in abundance during treatment. The MTX treatment is given at closer intervals than AP and MAP treatments, so there is always some drug at the tumor site during MTX treatment. Therefore, the changes in cell populations and cytokine concentrations over time for MTX regimen are smoother and do not fluctuate as much as in the other two treatments, even though the dynamics of MTX regimen follow the same trend as them. 

Since neoadjuvant chemotherapy tries to reduce the boundaries of the tumor before surgery, we can choose the desired size of tumor for surgery and run our optimization framework to find the optimal dosages of chemotherapy drugs for the tumor to reach this size at a specific time. Osteosarcoma sizes vary greatly among patients at first diagnosis, and large tumors cannot reduce to the same size as small tumors after neoadjuvant treatment without exceeding the safe dosages of chemotherapy drugs. Thus, we chose different desired cancer cell populations to optimize for depending on the size of tumor at first diagnosis. With MAP being the preferable treatment for osteosarcoma, as mentioned in the previous section, here we present the optimal dosages of the MAP regimen for large and small tumors in each cluster; the desired cancer cell population is 2.916 × 10 9 for large tumors and 1.36 × 10 8 for small tumors, which is equivalent to about 5 cm per dimension (length, width, depth) for large tumors and 1.8 cm per dimension for small tumors. We used 20,000 mg/m 2 of methotrexate, 45 mg/m 2 of doxorubicin, and 75 mg/m 2 of cisplatin per infusion day as maximum potential dosage, or equivalently 40,000 mg/m 2 of methotrexate, 90 mg/m 2 of doxorubicin, and 150 mg/m 2 of cisplatin per 35-day cycle.

The optimal dosages for large tumors are given in Table 5 . Large tumors were taken to be the steady state values of each cluster. As cluster 2 had the highest cancer cell population at the steady state, it had the highest optimal dosage for each drug of the MAP treatment among all clusters. Interestingly, clusters 1 and 3 had the same cancer cell population before treatment, but cluster 1 needed a higher dosage to achieve the same cancer cell population after the treatment as cluster 3. This was due to the fact that cluster 1's cancer cells grew faster during treatment, so the same dosage of drugs resulted in a higher cancer cell population in cluster 1 than in cluster 3 after treatment, as seen in Section 3.1. Thus, it is important to take the tumor growth rate of the patient into account while finding the optimal dosage of chemotherapy. Figure 8A shows that cancer cells in cluster 1 also grow fast after treatment, so it would be ideal to perform surgery quickly after neoadjuvant therapy. If it is impossible to start surgery promptly, we can choose a later time point to optimize for, so that at the time of surgery we still have the desired tumor size for resection. For example, if we cannot perform surgery until a month after chemotherapy, instead of using the cancer cell population at day 80 for optimization, which is 3 days after the second cycle of chemotherapy, we can use the cancer cell population at day 107 for optimization to find the optimal dosage, which is 30 days after the second cycle of chemotherapy. Then with the estimated optimal dosage, we will have the desired cancer cell population at day 107, which is the time of surgery.

The optimal dosages for a small tumor in each cluster are given in Table 6 . Small tumors were chosen with the same method as described in Section 3.4. The cancer cell populations in these small tumors were not much larger than the desired cancer cell population after treatment, so in all clusters the optimal dosages for small tumors are much smaller than the optimal dosages for large tumors. Cluster 3 particularly, whose cancer cell population before treatment was very close to the desired cancer cell population, has very small optimal dosages. In many cases, even though the tumor is small enough for resection, neoadjuvant chemotherapy is still given to remove any potential metastases that are too small to be yet detected. Another reason for neoadjuvant chemotherapy in small tumors is to allow evaluation of the tumor response [39] . Figure 8B suggests that although chemotherapy does not reduce cancer cell populations significantly as the cancer cell populations are already small to begin with, it helps prevent the cancer cell populations from growing bigger before surgery. Therefore, chemotherapy can also be used to control the growth ofa tumor while the patient waits for surgery. Overall, with our optimization framework, we can find an optimal chemotherapy dosage to obtain the desired cancer cell population on the day of surgery. Our results show that it is important to consider each individual patient's cancer growth rate while computing the optimal dosage, as patients with faster growth will need higher doses.

Cancer is a complex disease that consists of various components, such as immune cells, lymphatic vessels, and tumor cells [134] . Traditional in vivo and in vitro studies often explore one component of cancer at a time; thus, each study alone does not supply sufficient knowledge to understand cancer in all its heterogeneity, even though these studies do provide important insights about the individual cancer components they examine [135] . Mathematical modeling has also been utilized to study the components of cancer, especially the interactions between immune cells and cancer cells. However, the majority of the models only included one or two immune cells [136] [137] [138] [139] [140] [141] [142] [143] . Only a small subset of them have investigated the effects of multiple immune cells on cancer growth [52] [53] [54] . In addition, none of these studies applied their models to osteosarcoma. To the best of our knowledge, our previous work [58] was the first to study the interactions between cancer cells and the immune system in osteosarcoma.

With the increasing availability of gene expression data for many cancer types and the growing accuracy of tumor deconvolution methods, utilizing a deconvolution method on the gene expression data of a tumor to study the various components of the tumor microenvironment has become a more and more attractive option. Many recent studies applied the currently best performing deconvolution method, CIBERSORTx, to study the dynamics of cancer growth or to explore the relationships between clinical information and immune infiltrates [33, [144] [145] [146] . In this study, we extended our previous work [58] , which used CIBERSORTx to obtain immune abundances of osteosarcoma tumors and studied the tumor growth while considering its interactions with immune cells, to investigate the effects of chemotherapy on the osteosarcoma microenvironment.

Our results indicate that besides reducing the number of cancer cells, chemotherapy induces specific behaviors in certain immune cells and cytokines by causing necrosis of cancer cells. In particular, the populations of HMGB1 and dendritic cells increase when chemotherapy drugs are administered and decrease when these drugs are not given. In addition, helper T cells, cytotoxic cells, and IFN-γ generally increase in quantity during treatment, which aligns with the findings from [35, 36, 87, 90, [114] [115] [116] [117] [118] [119] [120] [121] [122] [123] [124] 126] . Meanwhile, cells and cytokines that are not affected by this necrotic cell death decrease in abundance due to being killed by chemotherapy drugs.

We note that it is good to start chemotherapy early, unless the tumor is close to its steady state, as tumors small or medium in size will grow more while the patient waits for treatment. It is especially important to start chemotherapy promptly for tumors that grow fast, such as those in cluster 1. Interestingly, we also noticed that even with the same initial cancer cell populations and the same dosages, the cancer cell population after treatment was higher in cluster 1 than in cluster 3 (cluster 3 has slower cancer growth rate than cluster 1). All of these observations suggest that it is necessary to take the unique growth rate of the tumor into consideration when choosing the dosage and treatment start time for the patient, thereby emphasizing the importance of personalized medicine.

In this study, we introduced a simple optimization framework to find the appropriate drug dosages to achieve a desired cancer cell population on a chosen day, such as the day of surgery. The results of our optimization also agree with the above observation that the individual's cancer growth rate is essential for calculating optimal chemotherapy dosages. Since high doses of chemotherapy are known to have high toxicity and to induce many serious health problems [147, 148] , it could be useful to use a mathematical model such as ours to estimate the appropriate dose rather than to give the standard dose for all tumor sizes, especially when small tumors are likely to need far smaller doses than the standard ones. Moreover, our model divides patients into groups based on their immune compositions, and thus can estimate their cancer growth more accurately than having one model for all patients, resulting in a more customized dosage recommendation for each patient.

Finding the right parameter values is a big challenge in the mathematical modeling of cancers. While it would be ideal to acquire parameters by performing in vivo and in vitro experiments, these experiments are often expensive and time-consuming. Therefore, numerous mathematical models approximate biologically reasonable values for parameters, or make assumptions about the relationships between parameters to derive their values, or vary the parameter values within certain ranges to study their effects on the outcomes. Here, we used chemotherapy-related parameters from a study that fitted these values to experimental data [84] , so our treatment parameters should be more accurate than those chosen based on biological rationality or derived from assumptions. For the two parameters that we had to assume appropriate values for, we studied their impacts on the results through sensitivity analysis and by varying them, and found that different values of these parameters result in fairly similar outputs with our model.

There are still some factors that our model does not account for. For instance, there are multiple levels of cancer cell sensitivity to chemotherapy, which means two different patients can both be resistant to a chemotherapy drug, but one patient might be more sensitive than the other. Thus, the drug efficacy coefficient of doxorubicin/cisplatinresistant cells used in our model does not represent all doxorubicin/cisplatin-resistant drug efficacies, as these parameters vary based on the levels of resistance of the cells. Despite that, our model is still useful for dose recommendations or for physicians to take into consideration while choosing between treatment options. Based on our model, a physician can monitor the tumor reduction throughout the treatment and adjust parameters such as drug efficacy coefficients according to how the tumor responds to treatment. However, it is important to note that although the results of our model align with some experimental observations in the literature, in vivo experiments should be performed to further validate our model before it can be utilized in clinic.

The chemotherapy-induced death rates of immune cells in our model were all approximated to be 0.6, so technically we could use one parameter to represent all of them. However, since chemotherapy targets cells with faster metabolic rates more successfully [89] , it is reasonable to expect that the death rates via chemotherapy differ between types of immune cells. Therefore, by keeping these death rates as separate parameters, our model can be easily improved by updating the rates of chemotherapy-induced death for immune cells in proportion to their metabolic rates. Other ways to improve upon this work include adding other chemotherapy drugs, such as ifosfamide, which is also a commonly used drug for osteosarcoma [39, 40] ; using partial differential equations to take into account the spatial distribution of the tumors as well; or extending it to different treatment options, such as radiotherapy and immunotherapy. To use our model for another type of treatment, one has to replace chemotherapy with the treatment of interest and include the interactions between that treatment and the tumor microenvironment.

In this study, we developed a mathematical model for the interaction network between the most common chemotherapy drugs and the key components of osteosarcoma microenvironment. We found that during chemotherapy, dendritic cells and HMGB1 increase in population when drugs are given and decrease in population while the patient waits for the next dose of drugs. Helper T cells, cytotoxic cells, and IFN-γ increase in abundance overall. Other cells and cytokines of the microenvironment that do not succumb to necrotic cell death have reduced populations after the treatment. Overall, the MAP regimen is effective at minimizing the number of cancer cells, and works better than methotrexate alone or a combination of doxorubicin and cisplatin. Our study also suggests the importance of considering the unique growth rate of the tumor when deciding on the dosage and the treatment start time for a patient, because fast growing tumors require higher dosages and earlier starts to treatment than slow growing tumors, as shown in our results. Data Availability Statement: Publicly available TARGET data were downloaded from UCSC Xena web portal: https://xenabrowser.net/datapages/ (accessed on 7 July 2021). Python scripts for computations and plotting the dynamical results are available here: https://github.com/ShahriyariLab/ Investigating-optimal-chemotherapy-options-for-osteosarcoma-patients-through-a-data-driven-mathemati (accessed on 7 July 2021).

The authors declare no conflict of interest.

The following abbreviations are used in this manuscript: 

(1 − e −β 1 A 1 ) + K M N (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 )

(1 − e −β 1 A 1 ) + K T h (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 )

(1 − e −β 1 A 1 ) + K T c (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 ) + K D N (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 ) + K D (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 ) + K C (1 − e −β 2 A 2 ) 

We introduce the following new non-dimensional variables:

Non-dimensional system (only equations with changes):

(1 − e −β 1 A 1 ) + K M N (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 )

(1 − e −β 1 A 1 )

(1 − e −β 1 A 1 ) + K T c (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 ) + K D N (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 ) + K D (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 ) + K C (1 − e −β 2 A 2 )

(1 − e −β 1 A 1 ) + K C (1 − e −β 2 A 2 )

where the new non-dimensional parameters are:

Key Statistics for Osteosarcoma

The Epidemiology of Osteosarcoma

Advances in limb salvage treatment of osteosarcoma

PDQ Pediatric Treatment Editorial Board. Osteosarcoma and Malignant Fibrous Histiocytoma of Bone Treatment (PDQ ® ): Patient Version; National Cancer Institute

Current Treatment Considerations for Osteosarcoma Metastatic at Presentation

Mechanisms of Resistance to Conventional Therapies for Osteosarcoma

A meta-analysis of randomized control trials of surgical methods with osteosarcoma outcomes

The Addition of Muramyl Tripeptide to Chemotherapy Improves Overall Survival-A Report From the Children's Oncology Group

Nivolumab in children and young adults with relapsed or refractory solid tumours or lymphoma (ADVL1412): A multicentre, open-label, single-arm, phase 1-2 trial

The Role of Radiotherapy in Oseosarcoma

Sternal Mass with Respiratory Compromise in a 10-year-old Child

Recent Insights into Therapy Resistance in Osteosarcoma

The tumor microenvironment in the response to immune checkpoint blockade therapies

Relevance of Immune Infiltration and clinical outcomes in pancreatic ductal adenocarcinoma subtypes

Tumoral immune-infiltrate (IF), PD-L1 expression and role of CD8/TIA-1 lymphocytes in localized osteosarcoma patients treated within protocol ISG-OS1

Tumor microenvironment and therapeutic response

Radiation fosters dose-dependent and chemotherapy-induced immunogenic cell death

Application of hyperthermia in addition to ionizing irradiation fosters necrotic cell death and HMGB1 release of colorectal tumor cells

The interaction between HMGB1 and TLR4 dictates the outcome of anticancer chemotherapy and radiotherapy

HMGB1-induced autophagy promotes chemotherapy resistance in leukemia cells

HMGB1 enhances immune suppression by facilitating the differentiation and suppressive activity of myeloid-derived suppressor cells

HMGB1: Guiding immunity from within

Emerging roles for HMGB1 protein in immunity, inflammation, and cancer

HMGB1: Endogenous danger signaling

Current and Emerging Targets in Immunotherapy for Osteosarcoma

Immunogenic Cell Death in Cancer Therapy

Innate immune cells: A potential and promising cell population for treating osteosarcoma

The role of interferons in the treatment of osteosarcoma

Regulatory T Cells and Cancer: A Two-Sided Story

Translational biology of osteosarcoma

Prognostic significance of HLA class I expression in osteosarcoma defined by anti-pan HLA class I monoclonal antibody

Immune Landscape of the Tumor Microenvironment Identifies Prognostic Gene Signature CD4/CD68/CSF1R in Osteosarcoma

Immune classification of osteosarcoma

Construction of immune-related gene pairs signature to predict the overall survival of osteosarcoma patients

Cisplatin-induced antitumor immunomodulation: A review of preclinical and clinical evidence

Tumor eradication by cisplatin is sustained by CD80/86-mediated costimulation of CD8+ T cells

Cisplatin alters antitumor immunity and synergizes with PD-1/PD-L1 inhibition in head and neck squamous cell carcinoma

Cisplatin augments cytotoxic T-lymphocyte-mediated antitumor immunity in poorly immunogenic murine lung cancer

Bone sarcomas: ESMO-PaedCan-EURACAN Clinical Practice Guidelines for diagnosis, treatment and follow-up

Chemotherapy and Other Drugs for Osteosarcoma

Symmetric vs. asymmetric stem cell divisions: An adaptation against cancer? PLoS ONE 2013, 8, e76195

The role of the bi-compartmental stem cell niche in delaying cancer

Modeling dynamics of mutants in heterogeneous stem cell niche

The role of backward cell migration in two-hit mutants' production in the stem cell niche

Mathematical models of cancer: When to predict novel therapies, and when not to

Hybrid modeling frameworks of tumor development and treatment

Cellular automaton models of tumor development: A critical review

Nonlinear modelling of cancer: Bridging the gap between cells and tumours

Cell dynamics in tumour environment after treatments

Implications of immune-mediated metastatic growth on metastatic dormancy, blow-up, early detection, and treatment

A stochastic model for cancer metastasis: Branching stochastic process with settlement

Mathematical modeling of regulatory T cell effects on renal cell carcinoma treatment

A model for effects of adaptive immunity on tumor response to chemotherapy and chemoimmunotherapy

A mathematical model of tumor-immune interactions

Mathematical modelling of the role of Endo180 network in the development of metastatic bone disease in prostate cancer

Polylactic acid/sodium alginate/hydroxyapatite composite scaffolds with trabecular tissue morphology designed by a bone remodeling model using 3D printing

A parameterised mathematical model to elucidate osteoblast cell growth in a phosphate-glass microcarrier culture

Data-Driven Mathematical Model of Osteosarcoma

Macrophages inhibit human osteosarcoma cell growth after activation with the bacterial cell wall derivative liposomal muramyl tripeptide in combination with interferon-γ

Tumor-associated macrophages in osteosarcoma: From mechanisms to therapy

TAMeless traitors: macrophages in cancer progression and metastasis

Immune Environment and Osteosarcoma

Distinct role of macrophages in different tumor microenvironments

The role of interleukin-2 during homeostasis and activation of the immune system

Plasticity of γδ T cells: Impact on the anti-tumor response

The contribution of immune infiltrates and the local microenvironment in the pathogenesis of osteosarcoma

The Osteosarcoma Microenvironment: A Complex However

TGF-β Signaling in Bone Remodeling and Osteosarcoma Progression

Guardian of T cell function

IL-10: The master regulator of immunity to infection

The two faces of IL-6 in the tumor microenvironment

Mesenchymal stem cells in the osteosarcoma microenvironment: Their biological properties, influence on tumor growth, and therapeutic implications

Emerging trends in immunotherapy for pediatric sarcomas

HMGB1 is an endogenous immune adjuvant released by necrotic cells

Necrosis of osteosarcoma cells induces the production and release of high-mobility group box 1 protein

HMGB1 in cancer: Good, bad, or both? Clin

Dendritic cells in the cancer microenvironment

How do regulatory t cells work?

T-cell-based immunotherapy for osteosarcoma: Challenges and opportunities

Interleukin 12 signaling in T helper type 1 (Th1) cells involves tyrosine phosphorylation of signal transducer and activator of transcription (Stat) 3 and Stat4

IL-12 produced by dendritic cells augments CD8+ T cell activation through the production of the chemokines CCL1 and CCL17

IL-12 stimulates CTLs to secrete exosomes capable of activating bystander CD8+ T cells

A mechanistic, predictive model of dose-response curves for cell cycle phase-specific and -nonspecific drugs

High-dose methotrexate with leucovorin rescue for severe COVID-19: An immune stabilization strategy for SARS-CoV-2 induced 'PANIC'attack

Doxorubicin: An update on anticancer molecular action, toxicity and novel drug delivery systems

Unexpected therapeutic effects of cisplatin

Immune response during therapy with cisplatin or radiation for human papillomavirus-related head and neck cancer

Mathematical model creation for cancer chemo-immunotherapy

Pretreatment with cisplatin enhances E7-specific CD8+ T-cell-mediated antitumor immunity induced by DNA vaccination

Cisplatininduced immune modulation in ovarian cancer mouse models with distinct inflammation profiles

A review of digital cytometry methods: Estimating the relative abundance of cell types in a bulk of cells

Primary tumor volume measurements in Ewing sarcoma: MRI inter-and intraobserver variability and comparison with FDG-PET

Size matters for sarcomas

Natural bone tissue and its biomimetic

Osteoblasts in bone tissue engineering

Simplified calculation of body-surface area

Determination of human body volume from height and weight

PubChem Compound Summary for CID 126941

PubChem Compound Summary for CID 443939, Doxorubicin Hydrochloride

PubChem Compound Summary for CID 5702198

Mixed immunotherapy and chemotherapy of tumors: Modeling, applications and biological interpretations

The Chemotherapy Source Book

Drugs & Diseases, Doxorubicin (Rx)

Cisplatin DrugBank Accession Number DB00515

Methotrexate DrugBank Accession Number DB00563

0: Fundamental algorithms for scientific computing in Python

Sensitivity analysis approaches applied to systems biology models

Likelihood approximation by numerical integration on sparse grids

Numerical integration using sparse grids

Comparison of MAPIE versus MAP in patients with a poor response to preoperative chemotherapy for newly diagnosed high-grade osteosarcoma (EURAMOS-1): An open-label, international, randomised controlled trial

Osteosarcoma MAP (Methotrexate, DOXOrubicin, cISplatin)

Neoadjuvant chemotherapy combined with limb salvage surgery in patients with limb osteosarcoma of Enneking stage II: A retrospective study

Co-delivery of doxorubicin and interferon-γ by thermosensitive nanoparticles for cancer immunochemotherapy

Immunological aspects of cancer chemotherapy

Immunogenic chemotherapy with cyclophosphamide and doxorubicin against established murine carcinoma

Dendritic cells combined with doxorubicin induces immunogenic cell death and exhibits antitumor effects for osteosarcoma

Immunogenicity of anthracyclines: Moving towards more personalized medicine

Doxorubicin Directs the Accumulation of Interleukin-12-Induced IFNγ into Tumors for Enhancing STAT1-Dependent Antitumor Effect

Mechanisms and management of doxorubicin cardiotoxicity

Combining tumor microenvironment modulating nanoparticles with doxorubicin to enhance chemotherapeutic efficacy and boost antitumor immunity

Lapatinib and doxorubicin enhance the S tat1-dependent antitumor immune response

Cisplatin augments antitumor T-cell responses leading to a potent therapeutic effect in combination with PD-L1 blockade

The mechanism of action of methotrexate

Caspase-dependent immunogenicity of doxorubicin-induced tumor cell death

Doxorubicin induces specific immune functions and cytokine expression in peritoneal cells

Silencing of immune activation with methotrexate in patients with COVID-19

Anti-inflammatory mechanisms of methotrexate in rheumatoid arthritis

Randomised trial of two regimens of chemotherapy in operable osteosarcoma: A study of the European Osteosarcoma Intergroup

Bone Cancer Treatment Regimens

Treatment of osteosarcoma of the extremities with the T-10 protocol, with emphasis on the effects of preoperative chemotherapy with single-agent high-dose methotrexate: A Scandinavian Sarcoma Group study

The efficacy and safety comparison of first-line chemotherapeutic agents (high-dose methotrexate, doxorubicin, cisplatin, and ifosfamide) for osteosarcoma: A network meta-analysis

Methotrexate, doxorubicin, and cisplatinum regimen is still the preferred option for osteosarcoma chemotherapy: A meta-analysis and clinical observation

Tumour heterogeneity and resistance to cancer therapies

Cancer Systems Biology: A peek into the future of patient care?

A mathematical model of antibody-dependent cellular cytotoxicity (ADCC)

Mathematical model of tumor-immune surveillance

The re-polarisation of M2 and M1 macrophages and its role on cancer outcomes

A dynamical model of tumour immunotherapy

Mathematical modelling for combinations of immuno-oncology and anti-cancer therapies

On study of immune response to tumor cells in prey-predator system

A model of dendritic cell therapy for melanoma

A validated mathematical model of tumor growth including tumor-host interaction, cell-mediated immune response and chemotherapy

Data Driven Mathematical Model of FOLFIRI Treatment for Colon Cancer

Immune classification of clear cell renal cell carcinoma

Data driven mathematical model of colon cancer progression

Prevention of High-Dose Chemotherapy-Induced Cardiotoxicity in High-Risk Patients by Angiotensin-Converting Enzyme Inhibition

Prevention and treatment of organ toxicity during high-dose chemotherapy: An overview. Anti-Cancer Drugs