Mathematical modeling of BCG-based bladder cancer treatment using socio-demographics

Cancer is one of the most widespread diseases around the world with millions of new patients each year. Bladder cancer is one of the most prevalent types of cancer affecting all individuals alike with no obvious “prototypical patient”. The current standard treatment for BC follows a routine weekly Bacillus Calmette-Guérin (BCG) immunotherapy-based therapy protocol which is applied to all patients alike. The clinical outcomes associated with BCG treatment vary significantly among patients due to the biological and clinical complexity of the interaction between the immune system, treatments, and cancer cells. In this study, we take advantage of the patient’s socio-demographics to offer a personalized mathematical model that describes the clinical dynamics associated with BCG-based treatment. To this end, we adopt a well-established BCG treatment model and integrate a machine learning component to temporally adjust and reconfigure key parameters within the model thus promoting its personalization. Using real clinical data, we show that our personalized model favorably compares with the original one in predicting the number of cancer cells at the end of the treatment, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$14.8\%$$\end{document}14.8% improvement, on average.


Introduction
Cancer is one of the most widespread illnesses in the world, responsible for millions of death every year with increasing numbers over time [1].Bladder Cancer (BC) is the seventh most common cancer worldwide, associated with 400 thousand new cases and 150 thousand deaths every year as of 2018 [2] and 600 thousand yearly new cases worldwide with only 77% five-year survival rate as of 2022 1 .BC has many forms and clinical stages, mainly differing by the depth of the cancer cell population in the urothelium [3].In the scope of this study, we focus on the non-invasive (superficial) BC where the cancer cells do not spread beyond the inner layer of the bladder where the entire cancer cell population is located inside the urothelium and does not invade other tissues.The non-invasive BC is highly common with roughly four out of five of all BC cases being diagnosed at the non-invasive stage [4].In these cases, multiple treatment protocols exist including chemotherapy-based [5] and immunotherapybased [6] treatments.Currently, the immunotherapy treatment suggested by [7] that follows weekly injections of Bacillus Calmette-Gérin (BCG) seems to achieve the best clinical improvement over a broad spectrum of clinical states [8,9].Most notably, BCG-based immunotherapy treatment has proven to reduce both the recurrence and progression of BC [10].The BCG treatment protocol is defined by the amount of the injected dosage, the number of injections, and the schedule of the treatment [11].Any change in one or more of these configurations can have a drastic effect on the patient's clinical state.However, due to the complexity of the biological dynamics, it is challenging to predict this change in advance.
In order to derive a suitable treatment protocol for patients, either at the individual or group level, researchers and clinicians often leverage the power of mathematical models and simulation [12].Commonly, in silico experiments provide a cheap, quick, and humane solution for clinical investigation of treatment protocols, allowing one to better understand and capture the underlying pharmacokinetics and pharmacodynamics [13,14].These models and simulations typically rely on an ordinary differential equation (ODE) representation where each variable describes a different cell population size [15,16,17,18,19].Indeed, the modeling and simulation of BC treatment protocols using this approach have attracted much attention in the literature [20,21,22].Notably, [23] proposed a BCG-based treatment protocol for BC which assumed continuous BCG instillation with a logistic growth for cancer cells inside the bladder.Formally, the proposed model takes the form: where B(t), E(t), T i (t), and T u (t) represent the concentration of BCG in the bladder, effector cell population size, the population of cancer cell that has been infected with BCG size, and the population of cancer cell that is uninfected with BCG size, respectively.The model's parameters represent the following quantities: p 1 is the rate of BCG killed by effector cells; p 2 is the infection rate of uninfected cancer cells by BCG; p 3 is the rate of destruction of cancer cell infected by BCG by effector cells; p 4 is the immune response activation rate; p 5 is the rate of effector cells deactivation after binding with infected cancer cells.α is the growth rate of effector cell population; λ is the cancer's population growing rate; b is the amount of BCG injected to the bladder.Fig. (1) shows a schematic view of this model, including the different cell populations and the interactions between them.One promising avenue for improving BCG-based treatment protocols is changing the currently practiced "onesize-fits-all" approach with a more personalized one [24,25].Specifically, by taking social and behavioral factors into account, one is likely to obtain a more favorable treatment prediction model at the individual patient's level and lead to better clinical outcomes as treatment optimization models would be able to have an underlined more accurate outcome prediction model [26].A patient's social and behavioral characteristics can be typically extracted with minimal overhead (e.g., by a simple questioning or directly from the patient's electronic health record) as opposed to alternative information-gathering efforts such as additional clinical tests which are associated with substantial operational costs.Following this line of thought, in this work, we propose a novel mathematical model which significantly extends that of [23].Our novelty lies primarily in the integration of a machine learning component which is used to assess and adjust the model's parameters over individuals and over time.Using real-world data of N = 417 patients, we show that our model favorably compares to the existing models.
The rest of the paper is organized as follows: Section 2 formally presents the framework, followed by Section 3 which provides theoretical outcomes regarding the proposed model.Next, Section 4 outlines the results of using the proposed model on real-world clinical data.Finally, in section 5, we analyze and discuss the results as well as propose possible future work directions.

Mathematical Modeling
Our model consists of two interconnected modules: a BCG-based treatment module and a socio-demographic personalization module.First, we define the bio-clinical dynamics of BCG treatment for BC.Then, we formalize the socio-demographics that underline the BCG-treatment dynamics.Based on these two modules, we formulate a fitting procedure to set the parameters of an instance of the framework using historical data.Fig. 2 shows a schematic view of the mathematical modeling.

BCG-based treatment module
Our following mathematical formulation relies on extensive prior literature which proposed and analyzed several biological models to describe the biological process underlying the BCG-based immunotherapy treatment for BC with increasing levels of complexity, capturing biological and clinical properties with great levels of detail and, presumably, accuracy [27,28,29,30,31,20].These and similar models describe, in a mathematical manner, the change in several cell populations over time due to (spatio-)temporal interaction between these cell populations [32,33,34].Generally speaking, the main line of work for modeling BCG-based treatment for BC, which we also follow in this work, was proposed by [21].The authors used a system of ODEs that represents the cell population sizes of several cell types over time.In particular, they divide the cell population into three main groups: BCG-infected, cancer, and immune-related cells, and described their interaction.In addition, special attention was placed on the distinction between BCG-infected and non-BCG-infected cells.
Here, we extend the model proposed by [23] in three manners: First, we replace the continuous BCG injection which assumes BCG is injected at some rate at any point in time with a discrete one which assumes a set of points in time in which BCG is injected.The latter more closely describes how BCG administration is provided in practice [21]; Second, we consider the uninfected cancer cell elimination by immune system cells [35]; Third, we introduce a healthy cell population and its interactions with the other cell types during a BCG treatment [25].Importantly, following [36]'s work, we assume that the BCG interaction with healthy and cancer cells is different following the difference between the cells' surfaces as well as shape.In addition, in order to allow personalization within the model, we replace the scalar parameters with functions that depend on time and the socio-demographics of the patient and add a term to the immune cells population that is associated with the presence of immune cells based on the level of the patient's activeness.Hence, the model takes the form (and explained right after): In Eq. ( 5), dB(t) dt is the dynamic change of BCG in the bladder over time.It is affected by the following five terms.First, a quantity b of BCG has been instilled into the bladder every τ steps in time.As the instillation of the BCG is modeled by a shifted Dirac delta function δ(t-mτ ), m ∈ {0, . . ., N -1}, the m th dose raises B(t) by b units at time t = mr.Second, BCG is eliminated by the immune cells (APCs) at a rate p 1 (t).Third and Fourth, BCG penetrates into uninfected cancer and uninfected healthy cells and is removed from the volume of the bladder while converting these cells into BCG-infected cancer cells at rates p 2 (t) and p 8 (t), respectively.Finally, the BCG cell population naturally decays at a rate µ B .
In Eq. ( 6), dE(t) dt is the dynamic number of immune cells over time.It is affected by the following five terms.First, the immune cell population naturally decays at a rate µ E (t).Second, immune cells are recruited due to the detection of BCG-infected cancer and healthy cells at a rate α(t).Third, immune cells are recruited due to bacterial infection in the bladder at a rate p 4 (t).Fourth and fifth, immune cells are destroyed while eliminating BCG-infected cancer and regular cells at rates p 5 and p 6 , respectively.
In Eq. ( 7), dTi(t) dt is the dynamic number of BCG-infected cancer cells over time.It is affected by the following two terms.First, BCG-infected cancer cells generated from uninfected cancer cells that interacted with BCG at a rate p 2 (t).Second, BCG-infected cancer cells are eliminated by immune cells at a rate p 3 (t).
In Eq. ( 8), dTu(t) dt is the dynamic number of uninfected cancer cells over time.It is affected by the following three terms.First, uninfected cancer cells naturally grow at a rate λ.Second, uninfected cancer cells become BCG-infected and are eliminated by immune cells at a rate p 2 (t).Third, immune system cells eliminate uninfected cancer cells at a rate p 3 .
In Eq. ( 9), dHu(t) dt is the dynamic number of uninfected healthy cells over time.It is affected by the following two terms.First, healthy cells are generated to fulfill the volume of the bladder, H m , at a rate p 7 (t).Second, uninfected healthy cells become BCG-infected due to the presence of BCG at a rate p 8 (t).
In Eq. ( 10), dHi(t) dt is the dynamic number of BCG-infected healthy cells over time.It is affected by the following two terms.First, uninfected healthy cells generated from uninfected cancer cells that interacted with BCG at a rate p 8 (t).Second, BCG-infected healthy cells are eliminated by immune cells at a rate p 9 (t).
For the proposed model, the initial condition at the beginning of the BCG treatment takes the form: A theoretical analysis of the model, proving that it is well-posed and analyzing its equilibria states and their stability, is provided in Section 3.

Socio-economic parameters
Patients are categorized into one of 72 socio-demographic groups.These groups are constructed based on the Cartesian product of four discretized properties: age (19-25, 26-35, 36-45, 46-55, 56-65, and 66+), gender (male and female), smocking behavior (smoker and non-smoker), and weight (underweight, normal weight, and overweight) such that the groups are pairwise disjoint.We decided on this democratization as it is commonly used in clinical studies and would be utilized later in this study as well [37,38,39,40].That said, the model is agnostic to the democratization of these parameters.Formally, each patient is represented by a timed finite state machine [41] as follows: p := (a, g, s, w) where a is the age group, g is the gender group, s is the smoking group, and w is the weight group.Importantly, it is assumed that a patient's socio-demographic properties do not change over the treatment process.

Treatment model fitting
Any model is as accurate as its fitting procedure allows.Thus, one may need to fit the proposed model on historical clinical data, obtaining a good approximation of the model's parameters' values.For the proposed model and in the clinical context it occurs, the historical data commonly have several limitations that make it challenging in the best case and infeasible in the worst case to use.Specifically, the BCG treatment data generally focuses on the treatment's clinical outcome and rarely, if any, include the cell population sizes (E, T i , T u , H i , and H u ) during the time of the treatment.As such, only two sample points during the course of the treatment are typically availableone at the beginning of the treatment and another one at its end.Second, the socio-demographic data is not directly taken into consideration in the model.However, it is correlated with the model's parameters.Third, the amount of available data is relatively small (usually, several hundred samples).Following the first point, and since measuring the amount of cancer cells is both clinically challenging and expensive, we assume that the data takes the form of T i (0) + T u (0), T i (t f ) + T u (t f ), ζ where t f ∈ N is the time at the end of the treatment and ζ ∈ R x is a vector of the socio-demographic properties of the patient such that x ∈ N is the number of socio-demographic properties.
During the fitting procedure, one is required to work with only the beginning and end point of the model which themselves are only providing partial knowledge of the model's state (i.e., T u and T i without B, E, H u , H m ).Hence, traditional fitting procedures might obtain unrealistic courses between these two points as long as the model closely crosses near them.To tackle this challenge, we proposed a three-step fitting procedure where each step is responsible to improve the accuracy and robustness of the model.First, we divide the data into k cohorts following the k-fold cross-validation method [42].This step is common in machine learning practices and improves the procedure's robustness [42].Then, each cohort is further divided into train, test, and validate sets.Using only the train set, we perform the first step.Namely, we adopt the fitting procedure proposed in [43] which, given the model's initial condition, the parameter space, historical data, and a loss function d, we utilize the gradient descent (GD) method [44] to find the parameters that minimize d on a fixed and finite duration in time [t 0 , t f ] such that t 0 < t f .Notably, the GD is applied on the model's parameters space such that the gradient for each configuration of parameter values is numerically obtained using the five-point stencil numerical scheme [45].The result of this process is the model's parameters' values that result in the closest clinical outcome of the given train set, divided into the socio-economic groups presented in the data.This step is shown to provide efficient initial fitting of the model's parameters [43].Nonetheless, due to the small amount of data, it is probably not spread across the domain properly which highly limits the usefulness of this method [46].To this end, for the second step, we randomly generate a new sample such that the features' values are range between the minimum and maximum of the values in the train set and chosen in a uniform distribution.Once a sample is generated, we test if the model's prediction of the clinical outcome's error is less or higher than the one obtained from a k nearest neighbor algorithm, on average, for the optimal choice of k.If the sample does not fulfill this condition, it is removed.Otherwise, it is added to the train set.This process repeats itself until n ∈ N samples are added.Intuitively, this step allows to generate a "filling" of synthetic samples based on the first step.Hence, afterward, the first step is repeated for the next train with a synthetic set.Using paired one-tail T-test, we check if the latter provides statistically significant better results on the test set.If it does not, we repeat the previous step.Otherwise, the train, synthetic, and test sets are merged and the last step is taking place.Namely, at this point, using the synthetic data from the second step, we remedy the shortcoming of the first step.Nonetheless, as we used the k nearest neighbor algorithm in the second step which is known to poorly extrapolate [47,48].Further improvement of the fitting robustness can be accomplished by utilizing other machine learning models for better extrapolation over the socio-demographic parameter space.Thus, for the third step, we use the Tree-based Pipeline Optimization Tool (TPOT) [49] automatic machine learning framework to search a large number of machine learning pipelines.For each pipeline, we test its performance on the validation set.The machine learning pipeline with the best performance, given a pre-defined metric (M ), is chosen.This process is repeated for each i ∈ [1, . . ., k] fold such that the overall prediction is the average of each machine learning model achieved at each fold.Overall, Fig. 3 provides a schematic view of the proposed fitting method and summarizes its main steps.

Theoretical Analysis
Here, we theoretically analyze the proposed model (Eq.5-10).We start by proving that the model always has a unique solution.Then, we identify the model's equilibria and analyze their stability.

Solution existence and uniqueness
In order to show that the proposed model has a solution and it is unique, we utilize the Picard-Lindelöf theorem [50].Formally, the Picard-Lindelöf theorem states that if D ⊂ R × R n is a closed rectangle with (t 0 , y 0 ) ∈ D and f : D → R n is a function that is continuous in t and Lipschitz continuous in y; then there exists some ϵ > 0 such that the initial value problem: has a unique solution y(t) on the interval [t 0 −ϵ, t 0 +ϵ].Thus, for our case y(t) := (B(t), E(t), T i (t), T u (t), H u (t), H i (t)).
In order to use the Picard-Lindelöf theorem, we first need to show that Eqs.(5-10) is continuous in t and Lipschitz continuous in y.To this end, let us consider a finite duration in time [0, T ] such that T < ∞.Next, the interaction between the of the unknown solution, y, has terms of a linear form and of the form y i y j , the function f such that dy(t)/dt = f (t, y(t) is C 1 which implies that it also locally satisfies Lipschitz condition and continuous in t [51].Therefore, one can apply the Cauchy-Lipschitz theorem [52] which leads to the result of the existence and uniqueness of the solution to Eq. (5-10), on any finite interval [0, T ].Hence, we show that a solution exists.As such, we need to show that for any non-negative initial condition, the solution is non-negative.To this end, let us assume a non-negative initial condition (B(0) ≥ 0, E(0) ≥ 0, T i (0) ≥ 0, T u (0) ≥ 0, H i (0) ≥ 0, H u (0) ≥ 0) for the proposed model (Eq.(5-10)).Now, let us focus on the fourth equation.By dividing by T u (t), one obtains: Computing the integral for t, we obtain that: Since T u (0) ≥ 0, we obtain T u (t) ≥ 0 for any value of t.For the fifth equation, we have a second-order equation in H u (t).Namely  takes the form: where a(t) = p 7 (t) − p 8 (t)B(t) − p7(t) Hi(t)+Tu(t)+Ti(t) Hm , b(t) = −p 7 (t)/H m .It is easy to see that as long as 0 ≤ p(t) than H u (t) ≥ 0 for any value of t.Using the fact that H u (t) > 0, we use the same method utilized for the third equation and obtain that 0 ≤ H i (t) for any t.The first and second equations yield that 0 ≤ E(t) and 0 ≤ B(t) for any t since T i (t), H i (t), T u (t), and H u (t) are non-negative.Thus, we show that the proposed model's solution is non-negative for any value of t > 0 if the initial condition is non-negative.

Equlibria and stability analysis
In order to better understand the bio-mathematical properties of the proposed model, we computed the equilibria states of the proposed model and their stability properties.Recall that an equilibrium state is reached when the system does not change without outside intervention.As such, to compute the equilibria states of the systems, we set the left side of the equations in Eqs.(5)(6)(7)(8)(9)(10) to zero and solve for the vector [B(t), Following this, one obtains two equilibria states: The first one is not mathematically valid, as for different values of t, B(t) would have different values during 0 ≤ t ≤ (N − 1)τ .As such, this equilibrium is well-defined for t > (N − 1)τ which results in 0 which is the trivial equilibrium where no dynamic takes place.For the second case, it is also trivial in the sense that all cells are healthy.
In order to obtain the equilibria states' stability of the two equilibria states, we first compute the Jacobian matrix of Eq. (??), following Routh-Hurwitz stability criterion [53]: ) Now, following the Hartman-Grobman theorem [54], by setting each equilibrium state to J we obtain: and We compute the eigenvalues of Eqs.(17)(18), obtaining that both have at least one eigenvalue that equals zero due to the third line being full of zeros, respectively.Hence, both equilibria are unstable.
4 Empirical Analysis

Data acquisition and preprocessing
We obtained retrospective data from the Bnai-Zion Medical Center 2 (Israel) from 2008 and 2017 [55].The data was anonymously extracted from the hospital's records under the following restrictions: 1) patients are adults 2 We refer the interested reader to https://www.gov.il/he/departments/b-zion-health-center/govil-landing-page (In Hebrew) (> 18 years old); 2) Patients received the standard BCG-based treatment for their non-invasive bladder cancer; 3) Patients were admitted between 2008 and 2017, in which period where all patients obtained the same treatment protocol for a non-invasive BC.In total, 417 BC patients are included, representing the entire patient population satisfying these three conditions.For each patient, we extracted the size of the cancer tumor at the beginning and end of the treatment alongside the four socio-demographic properties outlined above (gender, age, smocking, and weight).Additionally, for each patient, we extract the amount of BCG injected (b), the number of BCG injections (N ), the duration between every two consecutive BCG injections (τ ), BCG's decaying rate (µ B ), and the number of cells in the bladder (H m ).The BCG-treatment-related parameters (b, N, τ ) are defined by the standard treatment protocol [7] to be 2.8 • 10 8 , 6, and 7 days.The amount of cells in the bladder is highly linear to the patient's age and weight [56,57].As such, we use the formula proposed in [58] to approximate this value.The above is summarized in Table 1.
In addition, as T i (0) + T u (0) and T i (t f ) + T u (t f ) that required by the fitting procedure are not measurable "as is" rather than the polyp's volume it measured, one is required to map between these two values.To this end, we used the average volume of bladder cancer cells [59] and the polyp's volume to approximate T i (0) + T u (0) and T i (t f ) + T u (t f ).

H m
The number of healthy cells in the bladder without cancer as a function of the patient's gender and weight [1] 1.84 • 10 9 Table 1: The model's parameter definitions and average values as adopted from [24,25].

Parameter fitting
In order to use the proposed parameter fitting procedure, one is required to define a fitting metric.Hence, to make sure the model is able to better predict the outcome of a treatment procedure given an initial condition, we used the relative absolute error metric between the model's prediction and the historical data.Formally, let us denote the cancer population the model's prediction, and historical using T m t f and T h t f respectively.The relative mean absolute error than takes the form RM AE : ) such that P is the population of patients.Now, following the proposed fitting procedure (see Section 2.3) and the acquired data, we fitted the model and obtained 28.51 ± 3.76 relative absolute error on the test set of k = 5 folds.Fig. 4 shows the distribution of the relative mean absolute error among the four socio-demographic features -smoking habits, weight group, age group, and gender.We use black boxes to indicate configurations that could not be assessed given the available data (i.e., no single test case was available).Notably, non-smokers with underweight have the worst fitting results compared to other configurations.This may be partially attributed to the small number of samples within these groups during the training phase -13 in total.In addition, we compared the performance of the model's fitting, divided into smokers and non-smokers groups, using a two-tailed T-test [60].We obtained that zero is not included in the confidence interval of the T-test, which indicates that the two groups are not statistically significantly different.In a similar manner, when comparing the model's fitting performance between male to female patients using a T-test, we obtain p = 0.061 which indicates no statistically significant difference.Moreover, when comparing the weight categories using an ANOVA test with post-hoc Tukey pairwise tests [61], we see that the normal weighted category obtains statistically worse results compared to the two other categories p < 0.05.Similarly, the 56-65 and 65+ age groups obtained statistically significantly worse results compared to the younger age groups with p < 0.05.

Evaluation
We evaluate the proposed model by comparing it to two baselines: First, we use a narrow version of the proposed model which does not use of the socio-demographics at all.Second, we use the original model proposed by [23] as provided by the authors.Fig. 5 shows the relative mean absolute error of the three models for 50 samples.We choose the samples such that each sample would have a unique combination of the socio-demographic values.
One can clearly see that the proposed model outperforms both other baselines with an average relative absolute error of 19.38 ± 5.27 compared to 26.04 ± 6.85 and 34.18 ± 7.02, respectively.

Discussion
In this study, we presented a mathematical model to describe the dynamics of BCG-based immunotherapy for BC which incorporates socio-demographic parameters, such as age, gender, smoking behavior, and weight, to capture the heterogeneity in the patient populations and provide a more accurate treatment outcome prediction.By categorizing patients into distinct, easily identifiable, socio-demographic groups, we acknowledge the diverse characteristics of BC patients and their potential influence on the treatment outcome [62,63,64].This approach allows one to tailor more personalized models and thus obtain a better understanding of how different patient  profiles may respond to BCG-based immunotherapy.To this end, we used a novel fitting procedure that uses machine learning methods with only partial and sparse historical data to fit the personalized model for each sociodemographic group.
Our results strongly suggest that the currently unexploited socio-demographics may encompass clinically relevant information which underlines the BCG-based immunotherapy treatment of non-invasive BC.Most notably, the proposed model which takes these socio-demographics into account outperforms the baseline model proposed in [23] by 14.8.% as well as the same model which does not consider socio-demographics by 8.14%, as shown in Fig. 5. Arguably, these improvements should be attributed to the more personalized, thus accurate, ODE model.We observe further support for this observation in Fig. 4 which shows that different socio-demographic groups are associated with different fitting performances.This outcome can be associated with two factors: First, the differences in the amount of available data used to fit the model.Second, the clinical processes and dynamics that change as a result of belonging to each socio-demographic group.The effect of the first factor is clearly revealed for the non-smoker and the underweight group as a low amount of fitting data results in the highest, on average, fitting error.In a complementary manner, the second factor can be seen in older patients having higher fitting errors, on average, compared to younger patients while being a large portion of the dataset.
It is important to note that our model and results have certain limitations.First, the proposed model does not take into consideration the geometrical configuration of the bladder and therefore the spatial dynamics of the BCG treatment which has shown to have a critical role in the treatment outcomes and optimal treatment protocol for patients [25,33,65].As such, one can introduce these extensions to obtain a more realistic model.Second, the model assumes that a patient's socio-demographic properties remain constant during the course of BCG treatment.While this assumption is accepted due to the short duration of the treatment (i.e., around two months), it may not fully capture potential changes in patient characteristics that could influence treatment response.Third, the proposed model is partially personalized.Potentially, by using more parameters, one can obtain better personalization and even further improve the model's prediction accuracy [66].Finally, our model focuses on the treatment protocol clinical outcome prediction rather than suggesting an optimal personalized treatment as other studies do [24].In future work, we plan to further investigate this direction.
Overall, the proposed model provides a low (to no) overhead method for clinicians while providing a statistically significant improvement in the prediction of treatment outcomes.Thus, it can be easily adopted by healthcare professionals and help provide better treatment, saving and improving many lives.In the same manner, our proposed fitting method can be utilized in a broad spectrum of models as it assumes as little as possible on the historical data and nothing on the model itself.

Figure 1 :
Figure 1: Cell population dynamics in the bladder -taken with permission from [23].BCG (B) stimulates effector cells (E) of the immune system via APC activation.In addition, BCG infects uninfected cancer cells (T u ) which recruit effector cells into the bladder.Infected cancer cells (T i ) are destroyed by effector cells.

Figure 2 :
Figure 2: A schematic view of the mathematical modeling.
, for H u (t) the fifth equation is a Riccati equation, and therefore the solution of the equations and without the synthetic samples Using pariwise one-tail T-test, check if the train + synthetic datset provides statistically better results than the latter set for the Use automatic ML to test multiple ML pipelines Pick the model with the best performing for the validation set Average the parameter values for the k splits of the data

Figure 3 :
Figure 3: A schematic view of the proposed fitting procedure.

Figure 4 :
Figure 4: The average over k = 5 folds relative absolute error of the model's fitting on the test set.The black boxes indicate that no patients belong to this combination of socio-demographic values.

Figure 5 :
Figure5: A comparison between the proposed model, the proposed model without socio-demographics, and the standard model[23].