Mathematical model of STAT signalling pathways in cancer development and optimal control approaches

In various diseases, the STAT family display various cellular controls over various challenges faced by the immune system and cell death programs. In this study, we investigate how an intracellular signalling network (STAT1, STAT3, Bcl-2 and BAX) regulates important cellular states, either anti-apoptosis or apoptosis of cancer cells. We adapt a mathematical framework to illustrate how the signalling network can generate a bi-stability condition so that it will induce either apoptosis or anti-apoptosis status of tumour cells. Then, we use this model to develop several anti-tumour strategies including IFN-β infusion. The roles of JAK-STATs signalling in regulation of the cell death program in cancer cells and tumour growth are poorly understood. The mathematical model unveils the structure and functions of the intracellular signalling and cellular outcomes of the anti-tumour drugs in the presence of IFN-β and JAK stimuli. We identify the best injection order of IFN-β and DDP among many possible combinations, which may suggest better infusion strategies of multiple anti-cancer agents at clinics. We finally use an optimal control theory in order to maximize anti-tumour efficacy and minimize administrative costs. In particular, we minimize tumour volume and maximize the apoptotic potential by minimizing the Bcl-2 concentration and maximizing the BAX level while minimizing total injection amount of both IFN-β and JAK2 inhibitors (DDP).


Introduction
Needs extensive grammar revision. Eg, 'most fatal killer' is poor English. Paragraph 1 -"… two different dichotomous states of cell death: (i) an apoptosis progression and (ii) an anti-apoptosis progression." This is confusing. What is meant by anti-apoptosis progression? Does this also lead to cell death? If not, then how can it be a state of cell death? Paragraph 1 -What is meant by "apoptosis-mediated therapy"? This seems to be poorly worded.
Paragraph 2 -"However, there was no study of apoptosis mechanism through JAK-STAT signaling pathway with optimal control theory " This sentence has a couple of issues. First, it is poorly worded. Second, it is oddly specific. Is the claim that no mathematical studies of JAK-STAT mediation of apoptosis have been published? Or that no mathematical studies of JAK-STAT mediation of apoptosis using in some manner optimal control have been published? If the latter is the case, then this is hardly a novel study. In my experience, most mathematical models of cancer response to therapy include some manner of treatment optimization, whether or not they use optimal control. Where do the equations (1) -(4) come from? That is, why these functional form choices? Have they been derived from first principles? Are they purely phenomenological? How are they justified?
In general, it would be helpful to explain what the terms on the right hand side are doing rather than get into listing every parameter. Also, please include in the main text, a table with model variables and parameters and all units.
What is the difference between capital S (see paragraph 3 of introduction) and lowercase s (Page 5, lines [8][9]? What do you mean by IFN activity? This does not have any obvious mathematical meaning. Do you mean concentration of IFN? What units? In equation (11) -does the proliferation become negative when the apoptosis condition is satisfied? Then why have a separate apoptosis term? If not, then why does proliferation rate depend on apoptosis threshold? Should it not have it's own threshold? Where are the biological data justifying this? Is equation (11) also non-dimensional? Page 6, Lines 23-24 -what is meant by alternating injection, ie, alternating with what? By the way, this notation is very confusing. Stats are denoted by S_i and drug is also S? We hope you are keeping well at this difficult and unusual time. We continue to value your support of the journal in these challenging circumstances. If Royal Society Open Science can assist you at all, please don't hesitate to let us know at the email address below.

Dear Dr Kim
On behalf of the Editors, we are pleased to inform you that your Manuscript RSOS-210594 "Mathematical model of STAT signaling pathways in cancer development and optimal control approaches" has been accepted for publication in Royal Society Open Science subject to minor revision in accordance with the referees' reports. Please find the referees' comments along with any feedback from the Editors below my signature.
We invite you to respond to the comments and revise your manuscript. Below the referees' and Editors' comments (where applicable) we provide additional requirements. Final acceptance of your manuscript is dependent on these requirements being met. We provide guidance below to help you prepare your revision.
Please submit your revised manuscript and required files (see below) no later than 7 days from today's (ie 05-Aug-2021) date. Note: the ScholarOne system will 'lock' if submission of the revision is attempted 7 or more days after the deadline. If you do not think you will be able to meet this deadline please contact the editorial office immediately.
Please note article processing charges apply to papers accepted for publication in Royal Society Open Science (https://royalsocietypublishing.org/rsos/charges). Charges will also apply to papers transferred to the journal from other Royal Society Publishing journals, as well as papers submitted as part of our collaboration with the Royal Society of Chemistry (https://royalsocietypublishing.org/rsos/chemistry). Fee waivers are available but must be requested when you submit your revision (https://royalsocietypublishing.org/rsos/waivers).
Thank you for submitting your manuscript to Royal Society Open Science and we look forward to receiving your revision. If you have any questions at all, please do not hesitate to get in touch. Kind  Your paper is to be recommended for publications. Several comments of the referees, however, are very important and useful to make it improved. Please take regards.
Reviewer comments to Author: Reviewer: 1 Comments to the Author(s) The authors integrated a tumor growth model with intracellular signaling model. The model generated a bi-stability condition for cell state -apoptosis or anti-apoptosis. They considered two types of intramuscular injections (IFN-beta and cisplatin) as a medication to modulate the cell signal, which in turn controls the tumor growth. An optimal control theory was used to deduce optimal schedule to inject to reduce tumor growth more than constant and alternating treatment. I believe the findings of this study is a useful contribution to mathematical oncology field. A more detailed explanation of mathematical analysis and parameter estimation will be beneficial to readers.
Major Issues: 1. Model parameters are mentioned as estimated in the table and a detailed information is presented as supplementary. A brief explanation of sensitivity analysis and parameter estimation process in the main text will be useful for readers. 2. Mathematical details of the equilibria and bifurcation analysis would be helpful for readers to follow this study 3. How authors derived optimal control objective function equation (15)? What kind of numerical solution technique was used to solve the equation? 4. Between strategy 1 and 2, which one is better in what condition? 5. Discussion of implication, merit, a couple of limitation, future outlook would be helpful.
Reviewer: 2 Comments to the Author(s) This paper addresses intracellular signal transduction systems, including the STAT family. In particular, the control of the important bistable structures of cancer's apoptotic and anti-apoptotic states is investigated. This paper is the first to study the mechanism of apoptosis via the JAK-STAT signaling pathway using optimal control theory. As a therapeutic application, this study aims to find the optimal injection scheme of combination therapy (IFN-β + DDP) by optimal control theory. The model predicted that an optimally controlled schedule of these two drugs could provide better antitumor effects at minimal cost. Both the methods and results are new and valuable, and I recommend accepting this article for publication. The Introduction is a bit confusing. It will help eg at the end of the first paragraph to highlight knowledge gaps that mathematical modeling can fill. What is the need or motivation for the model?

Introduction
Needs extensive grammar revision. Eg, 'most fatal killer' is poor English. Paragraph 1 -"… two different dichotomous states of cell death: (i) an apoptosis progression and (ii) an anti-apoptosis progression." This is confusing. What is meant by anti-apoptosis progression? Does this also lead to cell death? If not, then how can it be a state of cell death?
Paragraph 1 -What is meant by "apoptosis-mediated therapy"? This seems to be poorly worded.
Paragraph 2 -"However, there was no study of apoptosis mechanism through JAK-STAT signaling pathway with optimal control theory " This sentence has a couple of issues. First, it is poorly worded. Second, it is oddly specific. Is the claim that no mathematical studies of JAK-STAT mediation of apoptosis have been published? Or that no mathematical studies of JAK-STAT mediation of apoptosis using in some manner optimal control have been published? If the latter is the case, then this is hardly a novel study. In my experience, most mathematical models of cancer response to therapy include some manner of treatment optimization, whether or not they use optimal control. Where do the equations (1) -(4) come from? That is, why these functional form choices? Have they been derived from first principles? Are they purely phenomenological? How are they justified?
In general, it would be helpful to explain what the terms on the right hand side are doing rather than get into listing every parameter. Also, please include in the main text, a Your revised paper should include the changes requested by the referees and Editors of your manuscript. You should provide two versions of this manuscript and both versions must be provided in an editable format: one version identifying all the changes that have been made (for instance, in coloured highlight, in bold text, or tracked changes); a 'clean' version of the new manuscript that incorporates the changes made, but does not highlight them. This version will be used for typesetting. Please ensure that any equations included in the paper are editable text and not embedded images.
Please ensure that you include an acknowledgements' section before your reference list/bibliography. This should acknowledge anyone who assisted with your work, but does not qualify as an author per the guidelines at https://royalsociety.org/journals/ethicspolicies/openness/.
While not essential, it will speed up the preparation of your manuscript proof if you format your references/bibliography in Vancouver style (please see https://royalsociety.org/journals/authors/author-guidelines/#formatting). You should include DOIs for as many of the references as possible.
If you have been asked to revise the written English in your submission as a condition of publication, you must do so, and you are expected to provide evidence that you have received language editing support. The journal would prefer that you use a professional language editing service and provide a certificate of editing, but a signed letter from a colleague who is a native speaker of English is acceptable. Note the journal has arranged a number of discounts for authors using professional language editing services (https://royalsociety.org/journals/authors/benefits/language-editing/).

===PREPARING YOUR REVISION IN SCHOLARONE===
To revise your manuscript, log into https://mc.manuscriptcentral.com/rsos and enter your Author Centre -this may be accessed by clicking on "Author" in the dark toolbar at the top of the page (just below the journal name). You will find your manuscript listed under "Manuscripts with Decisions". Under "Actions", click on "Create a Revision".
Attach your point-by-point response to referees and Editors at Step 1 'View and respond to decision letter'. This document should be uploaded in an editable file type (.doc or .docx are preferred). This is essential.
Please ensure that you include a summary of your paper at Step 2 'Type, Title, & Abstract'. This should be no more than 100 words to explain to a non-scientific audience the key findings of your research. This will be included in a weekly highlights email circulated by the Royal Society press office to national UK, international, and scientific news outlets to promote your work.

At
Step 3 'File upload' you should include the following files: --Your revised manuscript in editable file format (.doc, .docx, or .tex preferred). You should upload two versions: 1) One version identifying all the changes that have been made (for instance, in coloured highlight, in bold text, or tracked changes); 2) A 'clean' version of the new manuscript that incorporates the changes made, but does not highlight them. --If you are requesting a discretionary waiver for the article processing charge, the waiver form must be included at this step.
--If you are providing image files for potential cover images, please upload these at this step, and inform the editorial office you have done so. You must hold the copyright to any image provided.
--A copy of your point-by-point response to referees and Editors. This will expedite the preparation of your proof.

At
Step 6 'Details & comments', you should review and respond to the queries on the electronic submission form. In particular, we would ask that you do the following: --Ensure that your data access statement meets the requirements at https://royalsociety.org/journals/authors/author-guidelines/#data. You should ensure that you cite the dataset in your reference list. If you have deposited data etc in the Dryad repository, please only include the 'For publication' link at this stage. You should remove the 'For review' link.
--If you are requesting an article processing charge waiver, you must select the relevant waiver option (if requesting a discretionary waiver, the form should have been uploaded at Step 3 'File upload' above).
--If you have uploaded ESM files, please ensure you follow the guidance at https://royalsociety.org/journals/authors/author-guidelines/#supplementary-material to include a suitable title and informative caption. An example of appropriate titling and captioning may be found at https://figshare.com/articles/Table_S2_from_Is_there_a_trade-off_between_peak_performance_and_performance_breadth_across_temperatures_for_aerobic_sc ope_in_teleost_fishes_/3843624.

At
Step 7 'Review & submit', you must view the PDF proof of the manuscript before you will be able to submit the revision. Note: if any parts of the electronic submission form have not been completed, these will be noted by red message boxes.

Decision letter (RSOS-210594.R1)
We hope you are keeping well at this difficult and unusual time. We continue to value your support of the journal in these challenging circumstances. If Royal Society Open Science can assist you at all, please don't hesitate to let us know at the email address below.
Dear Dr Kim, I am pleased to inform you that your manuscript entitled "Mathematical model of STAT signaling pathways in cancer development and optimal control approaches" is now accepted for publication in Royal Society Open Science.
If you have not already done so, please remember to make any data sets or code libraries 'live' prior to publication, and update any links as needed when you receive a proof to check -for instance, from a private 'for review' URL to a publicly accessible 'for publication' URL. It is good practice to also add data sets, code and other digital materials to your reference list.
You can expect to receive a proof of your article in the near future. Please contact the editorial office (openscience@royalsociety.org) and the production office (openscience_proofs@royalsociety.org) to let us know if you are likely to be away from e-mail contact --if you are going to be away, please nominate a co-author (if available) to manage the proofing process, and ensure they are copied into your email to the journal. Due to rapid publication and an extremely tight schedule, if comments are not received, your paper may experience a delay in publication.
Please see the Royal Society Publishing guidance on how you may share your accepted author manuscript at https://royalsociety.org/journals/ethics-policies/media-embargo/. After publication, some additional ways to effectively promote your article can also be found here https://royalsociety.org/blog/2020/07/promoting-your-latest-paper-and-tracking-yourresults/.
On behalf of the Editors of Royal Society Open Science, thank you for your support of the journal and we look forward to your continued contributions to Royal Society Open Science. The authors integrated a tumor growth model with intracellular signaling model. The model generated a bi-stability condition for cell state -apoptosis or anti-apoptosis. They considered two types of intramuscular injections (IFN-beta and cisplatin) as a medication to modulate the cell signal, which in turn controls the tumor growth. An optimal control theory was used to deduce optimal schedule to inject to reduce tumor growth more than constant and alternating treatment. I believe the findings of this study is a useful contribution to mathematical oncology field. A more detailed explanation of mathematical analysis and parameter estimation will be beneficial to readers. Major Issues: 1. Model parameters are mentioned as estimated in the table and a detailed information is presented as supplementary. A brief explanation of sensitivity analysis and parameter estimation process in the main text will be useful for readers.
(Response): Thank you for careful reading and suggestions. For more explanation of sensitivity analysis, we replaced "Sensitivity analysis for twenty eight parameters in the model (Eqs ( (7)- (14)) was provided in Supplementary Information (SI) File." with "In order to see how sensitive are concentrations of main variables (STAT1, STAT3, Bcl-2, BAX, IFN-β, JAK2, DDP, and tumor) to twenty six parameters in the model (Eqs (9)-(16)) at various time points, we have performed sensitivity analysis. A partial rank correlation coefficient determines whether an increase (or decrease) in the parameter value will either decrease or increase the tumor volume and concentrations of main variables at a given time. See Supplementary Information File for more details." [ page 11, 1st paragraph, lines 1-6, in the (marked) revised manuscript ] For parameter estimation, we replaced "See Table I for parameter values in Eqs. (7)- (14). Parameter estimation was provided in Supplementary Information File." with "See Table I for parameter values in Eqs. (9)-(16) in a dimensionless form. Parameter values in the dimensional form are listed in Table S1 in the Supplementary Information File. Since our mathematical model contains many known and unknown parameter values, we provided parameter estimation in the Supplementary Information File, which is a necessary step toward fundamental and deep understanding of the dynamical process of the mathematical model. Parameter values are calculated based on empirical data such as half-life or estimated by fitting to experimental observation based on the mathematical structure of our model." [ page 6, 5th paragraph, lines 4-5 -page 7, 1st paragraph, lines 1-5, in the (marked) revised manuscript ] We also added the following parameter estimation in Supplementary Information File: "λ 1 , λ 2 : By using the reference value of Bcl-2 (B * = 10nM ), and taking the signal source of the Bcl-2, f 3 = 5.78 × 10 −2 nM/h, we estimate the dimensionless value of signaling strength to the Bcl-2 complex to be Similarly, by taking the reference value of BAX (X * = 351µM ), and taking the signal source of BAX, f 4 = 2.0288µM/h, we estimate the dimensionless value of signaling source to be λ 3 , k 5 : We assume the autocatalytic rate (a 7 ) of Bcl-2 is same as its negative leading to the dimensionless parameter value k 5 = a 7 µ S 1 B * = 1.0. By assuming that the signaling rate (λ ST AT 3 × (1.38µg/mL)) from STAT3 is at least 1.2-fold larger than its negative contribution (µ S 1 B * ) and by taking λ ST AT 3 = 2.513 × 10 −4 µM · mL/(h · µg), we get the dimensionless parameter value k 1 , k 3 , k 7 : We take the dimensionless parameter value k 1 = a 1 µ S 1 S * 1 = 4.0 by assuming that the autocatalytic rate (a 1 ) of STAT1 is at least 4-fold larger than its negative contribution (µ S 1 S * 1 ) from its decay in the absence of the inhibitory pathway from the STAT3. In a similar fashion, we get the autocatalytic rate of STAT3 (k 3 ) and BAX (k 7 ), k 2 , k 4 , k 6 , k 8 : The Hill-type coefficients a 2 , a 5 , a 8 , a 11 (k 2 , k 4 , k 6 , k 8 : parameters with dimensionless form) are fixed to be unity. k 9 , k 10 , µ T : In the absence of apoptosis, the carrying capacity of cancer is T 0 . However, in the presence of apoptosis, the carrying capacity is 10 +S 2 1 . Our system assumes that most cancer cells die by the apoptosis program. We set k 9 = 1.0, k 10 = 10, and µ T = 0.1 and the high level of STAT1 is about 4.5. Then, the steady state value of tumor volume in the presence of apoptosis is smaller than equilibrium value in the absence of apoptosis, The list of dimensionless parameters and their values are given in Table  I  The dynamics of the intracellular signaling system is governed by ordinary differential equations: In order to explore the equilibrium dynamics, the equilibrium state of the system must be identified first. This is done by setting Eqs. (1)-(4) to zero and solving for the solutions S E 1 , S E 3 , B E , and X E as follows: We can find the equilibrium point of the system (5)-(8) by nonlinear solver MATLAB as far as all parameter values are fixed. For example, we can find the equilibrium when S is 0, 0.4, and 1; when The Jacobian matrix, J , is given by: Then, the characteristic polynomial is given by where I is the identity matrix. We get two real eigenvalues (Λ 1 = −µ B < 0 and Λ 2 = −µ X < 0) and remaining two real eigenvalues, is smaller than 1, the eigenvalues Λ 3,4 are negative and the equilibrium is stable. On contrary, if Then, the characteristic equation is When Therefore, when S = 0, the equilibrium point is stable. When S = 1, AB; A ≈ 3.0568, B ≈ 0.0895, and AB ≈ 0.2736 < 1. There are three equilibrium points when S = 0.4. We calculated AB for each case. Two steady We can calculate AB and corresponding eigenvalues for continuous spectrum of IFN-β. Fig. S1A,B shows AB and Λ 3 for various IFN-β concentrations, respectively. The gray region represents the unstable region and stable region is marked in white. The green curve in Fig. S1A 1 corresponds to the lower branch in Fig. 3D in the main text and the red curve in unstable region in Fig. S1A 2 corresponds to the middle branch in Fig. 3D in the main text, lastly the blue curve in stable region in Fig. S1A 3 corresponds to the upper branch in Fig. 3D in the main text. The corresponding eigenvalues in these three regions are shown in Fig. S1B." [ Supplementary Information File, pages 3-5, Section 'Mathematical Analysis (equilibria and stability)', in the (marked) revised manuscript ] In addition, in order to make it more clear, we replaced the plain bifurcation STAT1 (S1)  Table I. 3. How authors derived optimal control objective function equation (15)? What kind of numerical solution technique was used to solve the equation?
(Response): Thanks for your valuable comment. We applied optimal control theory to find the "best" treatment regime that minimizes the tumor volume under the certain level of tumor concentration and the amount of the drugs used. The objective function in Eq. (15) consists of three parts: reducing the tumor volume in a certain level, remaining in an apoptosis state ({(B, X) : B < th B , X > th X }), and minimizing the use of drugs (u S and u D ). We used quadratic forms to simplify analysis with the convexity properties which are common in control problems of biological models [45]. For the controls in the integrand, we added linear terms to regularize the amount of drug used. In general, linear controls are more meaningful biologically than quadratic, but it is more difficult to analyze the system mathematically. We added the following phrase in the main text: "We used quadratic forms to simplify analysis with the convexity properties which are common in control problems in biological models [45]. For the controls in the integrand, we added linear terms to regularize the amount of drug used. In general, linear controls are more meaningful biologically than quadratic forms, but it is more difficult to analyze the system mathematically. Weight for each control is provided by ..." [ page 8, 1st paragraph, lines 8-12, in the (marked) revised manuscript ] To obtain the numerical solutions of the control problems, we used Forward-Backward Sweep Method (FBSM) which is based on shooting methods to solve boundary value problems [34]. We added the following statement at the end of the paragraph: "To obtain the numerical solutions of the control problems, we used forward-backward sweep method which is based on shooting methods to solve boundary value problems [34]." [ page 8, 1st paragraph, last two lines, in the (marked) revised manuscript ]

Between strategy 1 and 2, which one is better in what condition?
(Response): We are not comparing strategy 1 with strategy 2. While we in strategy 1 scheme, we were comparing optimal control strategy with other typical injection strategies (constant and periodic injection), we were comparing various injection combination with optimal control in a very specific case with 3 injections of IFN-β and 3 injections of DDP. In the latter case, we identified the best strategy out of many possible combination. It is hard to compare which one is better between strategy 1 and strategy 2 since it is in a different category.
5. Discussion of implication, merit, a couple of limitation, future outlook would be helpful.
(Response): Thank you for suggestions. For discussion of implication and merit, we added the following sentence in Conclusion Section "The mathematical model in this study may provide a comprehensive understanding of the IFN-β/JAK-induced STAT signaling network and the associated optimal control method may suggest an optimal infusion strategy of anti-cancer drugs in clinical setting." [ page 19, 1st paragraph, lines 3-6, in the (marked) revised manuscript ] For discussion of a couple of limitation, future outlook, we rewrote and added a paragraph in Conclusion Section as follows: "Our study has three main limitations: (i) Delivery of anti-cancer drugs is a complex process including transport of the anti-cancer agent through tissue. In this work, we did not take into account spatial movement of these agents. Therefore, a more general framework such as partial differential equations (PDEs) instead of the ODE model in this work may better describe the spatial transport of drugs. However, development of an optimal control scheme in the PDE model for a larger multi-scale system including blood vessels is still a challenge and we plan to investigate the spatial aspect of drug transport.
(ii) Our optimal control problems have not considered a linear form of controls which may be more realistic than a quadratic form. In particular, te ts u S (t)dt and te ts u D (t)dt represent the total amount of IFN-β and DDP, respectively. Therefore, this type of control forms in a minimization problem can be clearly interpreted as drug toxicity or cost by introducing weights [34,45,46]. We plan to develop a linear form of controls in a feasible setting of tumor models in future work. (iii) Tumor microenvironment and signaling networks play a major role in cancer progression and invasion. We did not take into account various microenvironmental factors including immune cells (macrophages, neutrophils, T cells, Th cells, T regs, and NK cells), cytokines/chemokines, and inter-and intra-cellular molecules. We plan to study these critical factors in future work.
Our mathematical formulation in this work can contribute to development of a new theoretical approach for other cell killing mechanism such as necroptosis [8,15,47] and autophagy [40] in cancer by the optimal control of key regulators in a ODE/PDE model [48,24,32,2,11] or multi-scale hybrid model [22,21,26,27,25,28,31,30] where key cellular death programs within the cancer cells can be taken into account at individual cell level in the multi-scale system." [ page 19, 3rd paragraph -6th paragraph, in the (marked) revised manuscript ] Minor issues: 1. The model is developed in equation (1) to (4) and then non-dimensionalized in (7) to (10). Further, the model is extended incorporating equations (11) to (14) without discussing the dimensional relevance of the variables and parameters appeared in (11) to (14) with the parameters in equation (6). We also added the following paragraph in order to provide the connections: "See Table I  (Response): In general, if controls have different scales, those weight constants could be different in order to keep the balance of effects on the objective function. For the same reason, C1 (C3) was chosen four times larger than C2 (C4) taking into account the amount of IFN-beta an DDP used. Meanwhile, the weight parameters C1 and C3 in our optimal problems are not necessarily different becasue, as we mentioned above, the linear terms are used for regularization of the controls.
4. T* is mentioned in the table but not appeared in the text. (Response): Thank you for careful reading and unclear expression in the caption of Fig 5. We meant that the down-regulation of STAT1/BAX and up-regulation of STAT3/BAX. Mathematically, it is an increasing function of S in the STAT1 graph for instance but we meant the STAT1 level stayed below the threshold. Same principle was applied to others. We now replaced "decreased activities of STAT1/BAX and increased levels of STAT3/Bcl-2." with "down-regulation of STAT1/BAX (i.e., S 1 < S th 1 , X < X th , ∀S (0 ≤ S ≤ 1)) and up-regulation of STAT3/Bcl-2 (i.e.,   This paper addresses intracellular signal transduction systems, including the STAT family. In particular, the control of the important bistable structures of cancer's apoptotic and anti-apoptotic states is investigated. This paper is the first to study the mechanism of apoptosis via the JAK-STAT signaling pathway using optimal control theory. As a therapeutic application, this study aims to find the optimal injection scheme of combination therapy (IFN-β + DDP) by optimal control theory. The model predicted that an optimally controlled schedule of these two drugs could provide better antitumor effects at minimal cost.
Both the methods and results are new and valuable, and I recommend accepting this article for publication. (5)] Should f 2 (s) be f 2 (s, j) ?

REFEREE 3:
In general the background referencing is incomplete. Researchers other than Kim et al. have done considerable work in modeling the Bcl-2 mediated regulation of apoptosis in cancer cells.
(Response): Thank you for careful reading and suggestions. We now added the following phrase in Introduction Section For example, mathematical models of Bcl-2 signaling networks illustrated importance of molecular play including intrinsic Bcl-2 apoptosis pathway [35,5,36,10,13], bistability in apoptosis [3], interaction between p53 and Bcl-2 [9], VEGF-Bcl-2 in angiogenesis [17,18], and MOMP regulation in pattern recognition [50]. See reviews in [52,19,16,42] for systems-based approaches of Bcl-2 and cell-death program. [ page 2, last paragraph, line 4 -page3, 1st paragraph, lines 1-4, in the (marked) revised manuscript ] The Introduction is a bit confusing. It will help eg at the end of the first paragraph to highlight knowledge gaps that mathematical modeling can fill. What is the need or motivation for the model? (Response): Thank you for suggestions. We restructured the Introduction Section and added the following phrases: Despite previous studies of apoptotic signaling, fundamental mechanism of the JAK-IFNβ-mediated apoptosis processes is poorly understood. However, translational studies with experimental data [38] support the considerable benefits of apoptosis-based therapy [12]. Qualitative analysis may contribute to fundamental understanding of this complex system. In particular, a new approach may identify the key functions and regulation of both JAK and IFNβ-mediated STATs in the apoptosis pathways within cancer cells. how changes in Bcl-2 and BAX affect cancer progression, " with "(i) unexplored structure of the STAT-JAK2-Bcl-2-BAX signaling pathways, (ii) how changes in IFN-β, STATs, and JAK2 affect cancer progression," [ page 3, 2nd paragraph, lines 4-5, in the (marked) revised manuscript ] We also added the following phrases at the end of the Introduction Section: "We found that JAK2 and mutual antagonism between STAT1 and STAT3 play a major role in regulation of the apoptosis and anti-apoptotic status in cancer cells, thus tumor growth dynamics, and obtained the optimal injection strategies of both JAK2 inhibitors and IFN-β by minimizing costs and maximizing antitumor efficacy through an optimal control theory. " [ page 3, 2nd paragraph, lines 6-9, in the (marked) revised manuscript ] So many parameters have been Est (see Table 1) that one cannot really have confidence in them. Estimated from where? Are these fit to some data? Why are most of them 1? (Response): Thank you for careful reading. Some of parameters such as decay rates are calculated from biological observation such as half-life of molecules. Many of estimated parameters came from fitting the simulation results to experimental data, for instance responses (up-and down-regulation) of intracellular molecules in response to periodic injection of IFN-β. Some of parameters are 1 because the mathematical structure of the mathematical model in a dimensionless form can gives us nice mathematical advantages to represent the biological behaviors such as switching of states in a given text. We now added the following in Supplementary Information File as well: "λ 1 , λ 2 : By using the reference value of Bcl-2 (B * = 10nM ), and taking the signal source of the Bcl-2, f 3 = 5.78 × 10 −2 nM/h, we estimate the dimensionless value of signaling strength to the Bcl-2 complex to be λ 1 = f 3 µ S 1 B * = 0.2. Similarly, by taking the reference value of BAX (X * = 351µM ), and taking the signal source of BAX, f 4 = 2.0288µM/h, we estimate the dimensionless value of signaling source to be λ 2 = f 4 µ S 1 X * = 0.2. λ 3 , k 5 : We assume the autocatalytic rate (a 7 ) of Bcl-2 is same as its negative contribution (µ S 1 B * = (2.89×10 −2 h −1 )×(1.0×10 −2 µM ) = 2.89×10 −4 µM/h), leading to the dimensionless parameter value k 5 = a 7 µ S 1 B * = 1.0. By assuming that the signaling rate (λ ST AT 3 ×(1.38µg/mL)) from STAT3 is at least 1.2-fold larger than its negative contribution (µ S 1 B * ) and by taking λ ST AT 3 = 2.513 × 10 −4 µM · mL/(h · µg), we get the dimensionless parameter value k 1 , k 3 , k 7 : We take the dimensionless parameter value k 1 = a 1 µ S 1 S * 1 = 4.0 by assuming that the autocatalytic rate (a 1 ) of STAT1 is at least 4-fold larger than its negative contribution (µ S 1 S * 1 ) from its decay in the absence of the inhibitory pathway from the STAT3. In a similar fashion, we get the autocatalytic rate of STAT3 (k 3 ) and BAX (k 7 ), k 3 = a 4 µ S 1 S * 3 = 4.0 and k 7 = a 10 µ S 1 X * = 4.0.
k 2 , k 4 , k 6 , k 8 : The Hill-type coefficients a 2 , a 5 , a 8 , a 11 (k 2 , k 4 , k 6 , k 8 : parameters with dimensionless form) are fixed to be unity. k 9 , k 10 , µ T : In the absence of apoptosis, the carrying capacity of cancer is T 0 . However, in the presence of apoptosis, the carrying capacity is T 0 1 − µ T r 1− k 9 S 2 1 k 2 10 +S 2 1 .
Our system assumes that most cancer cells die by the apoptosis program. We set k 9 = 1.0, k 10 = 10, and µ T = 0.1 and the high level of STAT1 is about 4.5. Then, the steady state value of tumor volume in the presence of apoptosis is smaller than equilibrium value in the absence of apoptosis, The list of dimensionless parameters and their values are given in Table I  (Response): Thank you for nice suggestions. We now replaced "The mathematical model unveils the structure and functions of the intracellular signaling and cellular outcomes of the anti-tumor drugs in the presence of IFN-β and JAK stimuli. We finally use an optimal control theory in order to suggest optimal anti-tumor efficacy with minimal costs in response to drug administration." with "The roles of JAK-STATs signaling in regulation of the cell death program in cancer cells and tumor growth are poorly understood. The mathematical model unveils the structure and functions of the intracellular signaling and cellular outcomes of the anti-tumor drugs in the presence of IFN-β and JAK stimuli. We identify the best injection order of IFN-β and DDP among many possible combinations, which may suggest better infusion strategies of multiple anti-cancer agents at clinics. We finally use an optimal control theory in order to maximize anti-tumor efficacy and minimize administrative costs. In particular, we minimize tumor volume and maximize the apoptotic potential by minimizing the Bcl-2 concentration and maximizing the BAX level while minimizing total injection amount of both IFN-β and JAK2 inhibitors (DDP)." [ Abstract, lines 7-17, in the (marked) revised manuscript ] 2.
[Abstract] Line 11 what specifically is meant by immunity?
(Response) In order to make it clear, we replaced "immunity" with "various challenges faced by the immune system" [ Abstract, lines 1-2, in the (marked) revised manuscript ] 3.
(Response): Thank you for careful reading. We now replaced "cellular fate" with "cellular states" 5.
[Introduction] Paragraph 1 What is meant by apoptosis-mediated therapy? This seems to be poorly worded.
(Response): Thank you for pointing out this. More commonly used term is 'apoptosis-based therapy'. We now replaced "apoptosis-mediated therapy" with "apoptosis-based therapy [12]" [ page 2, 1st paragraph, line 22, in the (marked) revised manuscript ] in the edited phrases in 1st paragraph as following: "Despite previous studies of apoptotic signaling, fundamental mechanism of the JAK-IFNβ-mediated apoptosis processes is poorly understood. However, translational studies with experimental data [38] support the considerable benefits of apoptosis-based therapy [12]. Qualitative analysis may contribute to fundamental understanding of this complex system. In particular, a new approach may identify the key functions and regulation of both JAK and IFNβ-mediated STATs in the apoptosis pathways within cancer cells." [ page 2, 1st paragraph, lines 20-25, in the (marked) revised manuscript ] 6.
[Introduction] Paragraph 2 'However, there was no study of apoptosis mechanism through JAK-STAT signaling pathway with optimal control theory.' This sentence has a couple of issues. First, it is poorly worded. Second, it is oddly specific. Is the claim that no mathematical studies of JAK-STAT mediation of apoptosis have been published? Or that no mathematical studies of JAK-STAT mediation of apoptosis using in some manner optimal control have been published? If the latter is the case, then this is hardly a novel study. In my experience, most mathematical models of cancer response to therapy include some manner of treatment optimization, whether or not they use optimal control.
(Response): Thank you for pointing out this. It's indeed the latter case. We replaced "In particular, optimal control approaches are used to identify optimal schedule of anti-cancer drugs [2,20,43,44]. However, there was no study of apoptosis mechanism through JAK-STAT signaling pathway with optimal control theory." with "In particular, optimal control approaches are used to identify optimal schedule of anti-cancer drugs targeting stromal/immune cells and various signaling pathways [2,20,43,44]. Fundamental mechanism of the JAK-STAT-mediated cancer cell killing is still poorly understood. (Response): In mathematics, we typically use variable symbols with an overbar in oder to indicate they are variables with dimension (sometimes dimensionless variables depending on authors and styles). Sometimes, some researchers use a 'tilde' instead of an overbar. As we explained below, we are using many symbols for dimensional or dimensionless variables and parameters. Therefore, this was our choice to make it more clear rather than making it confused by introducing so many different symbols. If we introduce too many symbols without using bar, we would be running out of symbols and making it confusing.

[Methodology]
Where do the equations (1) -(4) come from? That is, why these functional form choices? Have they been derived from first principles? Are they purely phenomenological? How are they justified?
(Response): The equations (1) -(4) are based on mass balance equation. Since, we don't use partial differential equation, the simple mass balance concept 'input -output = accumulation' concept is used, i.e., dy dt = In -Out. Let's supposed that we have a system of ODEs with N variables (y i = y i (t), (i = 1, . . . , N )). In general, the mass balance of given intracellular variable y i = y i (t), (i = 1, . . . , N ) is used to derive governing equation where y = (y 1 , y 2 , . . . , y N ), the function f i (y) represents the source, g i (y) represents inhibition, and h i (y) represents outflux due to natural decay, i.e. h i (y) = µ i y i , where µ i is the decay rate. The source function f i (y) can be described below based on biological observations. A fractional form for the inhibition term in Eq (17) was chosen as the qualitative representation of negative feedbacks in this work. Specifically, we use the form for autocatalytic activity with the inhibition process of the intracellular variable y i by another intracellular variable y j (i = j), where ζ 1 , ζ 2 are constants, the parameter α i represent the inhibition strength along with amount of the variable y j via a function F (y i ) (µ i , ζ 1 , ζ 2 , α i ∈ R + , n ∈ Z + ). In the absence of source, this inhibition term with the decay term,−µ i y i provides the baseline concentration y * i ≈ ζ 1 µ i of the given molecule y i at a steady state when the inhibition strength F (y j ) is small or zero. (When f i (y) = 0, the baseline becomes y * i ≈ f i +ζ 1 µ i .) The relative balance between the source term and inhibition strength from y j essentially determines the concentration of the molecule y i . Thus, by comparing the simulated y i level with experimental data in the presence and absence of the inhibitory molecule y j in the system, one can build a mathematical model in Eq (17) with the consistent, up-or down-regulated y i . Several studies [1,29,23,22] have shown that this fractional form for the negative feedbacks may reproduce analytic structure of genetic networks (positive and negative feedbacks) and qualitative dynamics such as bi-stability with experimental validation. Other forms of negative feedbacks (eg. one based on chemical reactions) have been used in the literature [49,25]. Essentially, this choice allows us to model the inhibition process between intracellular molecules, which was observed in experiments. Other forms are possible, of course. However, other choices typically involves countless chemical reactions and another choices of model design based on chemical reactions at some point when one wants to model the negative feedback of the given molecule by another, which generates same degrees of uncertainty, in other words every model design have their own assumptions. In order to make it clear, we now replaced "The mass balance.." with "In this work, we ignore any spatial effects on dynamics of given system. In general, the mass balance of given intracellular variable y i = y i (t), (i = 1, . . . , N ) is used to derive governing equation where y = (y 1 , y 2 , . . . , y N ), the function f i (y) represents the source, g i (y) represents inhibition, and h i (y) represents outflux due to natural decay, i.e. h i (y) = µ i y i , where µ i is the decay rate. The source function f i (y) can be described below based on biological observations. A fractional form for the inhibition term in Eq (17) was chosen as the qualitative representation of negative feedbacks in this work. Specifically, we use the form for autocatalytic activity with the inhibition process of the intracellular variable y i by another intracellular variable y j (i = j), where ζ 1 , ζ 2 are constants, the parameter α i represent the inhibition strength along with amount of the variable y j via a function F (y i ) (µ i , ζ 1 , ζ 2 , α i ∈ R + , n ∈ Z + ). In the absence of source, this inhibition term with the decay term,−µ i y i provides the baseline concentration y * i ≈ ζ 1 µ i of the given molecule y i at a steady state when the inhibition strength F (y j ) is small or zero. (When f i (y) = 0, the baseline becomes y * i ≈ f i +ζ 1 µ i .) The relative balance between the source term and inhibition strength from y j essentially determines the concentration of the molecule y i . Thus, by comparing the simulated y i level with experimental data in the presence and absence of the inhibitory molecule y j in the system, one can build a mathematical model in Eq (17) with the consistent, up-or down-regulated y i . Several studies [1,29,23,22] have shown that this fractional form for the negative feedbacks may reproduce analytic structure of genetic networks (positive and negative feedbacks) and qualitative dynamics such as bi-stability with experimental validation. Other forms of negative feedbacks (eg. one based on chemical reactions) have been used in the literature [49,25]. Then, the mass balance.  The autocatalytic activity parameters for the modules of STAT1, STAT3, Bcl-2, and BAX are denoted by a 1 , a 4 , a 7 , and a 10 , respectively. The parameters a 2 , a 5 , a 8 , and a 11 are the Hill-type inhibition saturation constants from the counter part of STAT1, STAT3, Bcl-2, and BAX, respectively. There are constant of JAK. u D will be set similarly to u S ." with "JAK2 in Eq (15) undergoes production at a rate J s , the DDP-mediated degradation of JAK by DDP (γ D ), and decay process at a rate µ J in the first, second, and third terms on RHS, respectively. The first and second terms on RHS in Eq (16) represent the time-dependent injection of DDP via a function u D (t) and decay process at a rate µ D , respectively." [ page 6, 4th paragraph, lines 7-10, in the (marked) revised manuscript ] Yes, 'IFN activity' means 'high concentration of IFNβ' with a typical unit (g/mL or µM ) (these two words are almost synonym in biology since it is sometimes hard to measure concentration of intracellular (or extracellular) molecules). In order to make it clear, we made the following changes: (Response): Thank you for careful reading and suggestions. We replaced "another one (th X ) in genetic classification, in other words, when the mathematical condition {(B, X) : B < th B , X > th X } is satisfied. The threshold values were set based on biological observation and dynamical system of Eq. (9)- (12)." with "another one (th X ), in other words, when the condition {(B, X) : B < th B , X > th X } is satisfied, as suggested in experiments and modeling works [14,4,33,3,51,39,24]. The threshold values were set based on biological observation [14,4,33,3,51,39] and dynamical system of Eq. (9)- (12)." [ page 5, 3rd paragraph, lines 4-7, in the (marked) revised manuscript ] 13.
[Methodology] In equation (11) does the proliferation become negative when the apoptosis condition is satisfied? Then why have a separate apoptosis term? If not, then why does proliferation rate depend on apoptosis threshold? Should it not have its own threshold? Where are the biological data justifying this? Is equation (11) also non-dimensional?
(Response): Yes, the equation (11) (now Eq (13) in the revised version) is non-dimensional. The nondimensionalization of these variables (tumor volume, IFN-β, DDP, JAK) was provided in Supplementary Information File. In equation (11) (now Eq (13) in the revised version), the proliferation does not become negative when the apoptosis condition is satisfied. Even when the proliferation term is zero in the first term on the right hand side (RHS), for example when r = 0, when the apoptosis condition is satisfied (i.e., when I {B<th B , X>th X } = 1), the equation becomes dT dt = − µ T T apoptosis . and the solution of the equation is T (t) = χ 0 e −µ T t with the initial tumor size k 9 S 2 1 k 2 10 +S 2 1 I {B<th B , X>th X } ≥ 0, ∀S 1 due to our assumption k 9 ≤ 1. In particular, we set k 9 = 1." [ page 6, 2nd paragraph, lines 8-10, in the (marked) revised manuscript ] Inhibition of tumor growth by STAT1 in TME was observed in literature including [6] as described in the paragraph. We have the separate apoptosis term since tumor cells can be killed under the apoptosis condition even when the STAT1 level is low. Switching of the STAT1-mediated inhibition of tumor growth is modeled by adapting the Hill type function k 9 S 2 1 k 2 10 +S 2 1 rather than using threshold value of STAT1. This Hill type function has similar effect as threshold values. Therefore, those two mechanisms generate combined tumor cell killing in our modeling framework.
(i) Delivery of anti-cancer drugs is a complex process including transport of the anti-cancer agent through tissue. In this work, we did not take into account spatial movement of these agents. Therefore, a more general framework such as partial differential equations (PDEs) instead of the ODE model in this work may better describe the spatial transport of drugs. However, development of an optimal control scheme in the PDE model for a larger multi-scale system including blood vessels is still a challenge and we plan to investigate the spatial aspect of drug transport. (ii) Our optimal control problems have not considered a linear form of controls which may be more realistic than a quadratic form. In particular, te ts u S (t)dt and te ts u D (t)dt represent the total amount of IFN-β and DDP, respectively. Therefore, this type of control forms in a minimization problem can be clearly interpreted as drug toxicity or cost by introducing weights [34,45,46]. We plan to develop a linear form of controls in a feasible setting of tumor models in future work. (iii) Tumor microenvironment and signaling networks play a major role in cancer progression and invasion. We did not take into account various microenvironmental factors including immune cells (macrophages, neutrophils, T cells, Th cells, T regs, and NK cells), cytokines/chemokines, and inter-and intra-cellular molecules. We plan to study these critical factors in future work.
Our mathematical formulation in this work can contribute to development of a new theoretical approach for other cell killing mechanism such as necroptosis [8,15,47] and autophagy [40] in cancer by the optimal control of key regulators in a ODE/PDE model [48,24,32,2,11] or multi-scale hybrid model [22,21,26,27,25,28,31,30] where key cellular death programs within the cancer cells can be taken into account at individual cell level in the multi-scale system." [ page 19, 3rd paragraph -6th paragraph, in the (marked) revised manuscript ]