Optimization of SOX2 Expression for Enhanced Glioblastoma Stem Cell Virotherapy

TL;DR

This study uses a mathematical model to show how adjusting SOX2 levels can improve Zika virus treatment for brain cancer stem cells.

Contribution

A novel ODE model reveals critical thresholds for SOX2 and viral replication to optimize Zika virotherapy for glioblastoma.

Findings

Critical thresholds for SOX2 and viral replication were identified to change equilibrium stability.

Transcritical and Hopf bifurcations were observed in population dynamics affecting therapeutic outcomes.

Optimal SOX2 expression and viral bursting size enhance Zika virotherapy efficacy against glioblastoma stem cells.

Abstract

The Zika virus has been shown to infect glioblastoma stem cells via the membrane receptor αvβ5, which is activated by the stem-specific transcription factor SOX2. Since the expression level of SOX2 is an important predictive marker for successful virotherapy, it is important to understand the fundamental mechanisms of the role of SOX2 in the dynamics of cancer stem cells and Zika viruses. In this paper, we develop a mathematical ODE model to investigate the effects of SOX2 expression levels on Zika virotherapy against glioblastoma stem cells. Our study aimed to identify the conditions under which SOX2 expression level, viral infection, and replication can reduce or eradicate the glioblastoma stem cells. Analytic work on the existence and stability conditions of equilibrium points with respect to the basic reproduction number are provided. Numerical results were in good agreement with analytic solutions. Our results show that critical threshold levels of both SOX2 and viral replication, which change the stability of equilibrium points through population dynamics such as transcritical and Hopf bifurcations, were observed. These critical thresholds provide the optimal conditions for SOX2 expression levels and viral bursting sizes to enhance therapeutic efficacy of Zika virotherapy against glioblastoma stem cells. This study provides critical insights into optimizing Zika virus-based treatment for glioblastoma by highlighting the essential role of SOX2 in viral infection and replication.

Figures

A](#F2),[B](#F2)). The equilibrium point $[eqn]$ becomes asymptotically stable when $[eqn]$ , while $[eqn]$ becomes asymptotically unstable and it satisfies the condition $[eqn]$ . The solution $[eqn]$ converges to $[eqn]$ . This result provides a good agreement between our analytic and numerical results.

Equations

Full text

Introduction

Glioblastoma multiforme (GBM), commonly known as glioblastoma, is one of the most frequent and malignant of all central nervous system (CNS) tumors [1–3]. This aggressive cancer currently leads to a life expectancy of 15 months once diagnosed, with fewer than 5% of patients surviving for five years after diagnosis [4–6]. Glioblastoma prognosis is extremely poor due to its resistance to treatment and its aggressive nature [7,8]. Patients can experience a diverse range of symptoms, such as increased intracranial pressure, headaches, focal or progressive neurologic deficits, and seizures [2]. The treatments used for malignant primary brain tumors are radiotherapy with concomitant chemotherapy using the alkylating agent temozolomide and surgical resection [2,8,9]. Despite available treatments to inhibit tumor growth, GBM remains incurable [3]. According to the World Health Organization (WHO), GBM is classified as a grade IV glioma, the highest grade [6,10,11].

Cancer stem cells are known for their pluripotent and self-renewing properties, contributing to the initiation and progression of tumors [7]. Glioblastoma stem cells (GSCs) play an important role in mediating therapeutic resistance and recurrence in glioblastoma [7,9]. If not eliminated during chemotherapy, GSCs tend to reinitiate tumor formation, causing tumor recurrence [12]. Therefore, the development of novel therapeutic strategies to specifically target GSC populations for elimination could have significant clinical impacts [7].

In oncolytic virotherapy (OVT), lytic viruses are engineered to specifically target cancer cells, leading to viral replication within cancer cells and host-cell lysis [13,14]. In contrast to conventional radiation therapy or chemotherapy, OVT minimizes damage to normal cells and tissues, offering a more targeted and tolerable therapeutic approach [13–15]. The first FDA-approved virotherapy, Talimogene laherparepvec (Amgen Inc., Thousand Oaks, CA, USA), utilizes an engineered herpes simplex virus to target and lyse melanoma cancer cells [16]. Various types of oncolytic viruses have been employed to treat GBM. Among these, the mutated herpes virus G47 has been clinically used for recurrent glioblastoma, with its usage approved in Japan [17]. Recently, Ling et al. reported their phase I trial data involving 41 GBM patients, demonstrating that their engineered herpes virus CAN-3110 could induce T cells and improve patient survival [18].

Zika virus (ZIKV) belongs to the flavivirus genus of RNA viruses and infects the central nervous system by targeting neural precursor cells, ultimately leading to cellular death [19]. ZIKV is known to cause neurodevelopmental disruptions and brain abnormality in fetuses, including microcephaly [19]. While ZIKV can replicate in adult brain tissue and target mature neurons [20], the symptoms following viral infection in adults generally appear to be minimal [21]. Due to the differentiation, proliferation, and cell death observed during ZIKV infection, coupled with the smaller number of negative effects seen in adults, researchers have explored the potential use of ZIKV as a virotherapeutic agent in various neural cancers, including glioma [22,23], neuroblastoma [24,25], and glioblastoma [26–29].

Several mathematical models have been developed to elucidate the intricate interactions between cancer cells and oncolytic viruses employing ordinary differential equations (ODEs) [30–33] and partial differential equations (PDEs) [19,34,35], incorporating time delays to address the lytic cycle of viruses [34–36]. Furthermore, studies have proposed the anti-viral and anti-tumoral effects of immune responses on OVT [37–40]. Our previous study was to investigate the role of natural killer (NK) cells on OVT where equilibrium points can be created or destroyed by activation of NK cells [37]. Recent experimental studies [27–29]show that ZIKV selectively infects GSCs via the membrane receptor $[eqn]$ , which is highly expressed by the transcription factor SOX2 in GSCs. This positive correlation between the levels of SOX2 expression and ZIKV infection suggests the potential utility of SOX2 as a therapeutic marker for predicting successful treatment outcomes. This also underscores the significance of monitoring SOX2 levels during ZIKV virotherapy to optimize therapeutic approaches. Although the complex interplays between cancer cells and oncolytic viruses have been investigated, no mathematical model has been developed to address the dynamic interaction between cancer stem cells and Zika viruses.

In this study, we developed a mathematical ODE model to understand the underlying mechanism of interaction among GSCs, infected GSCs, and ZIKV. We focused on examining the effect of SOX2 expression levels and various viral bursting sizes on the outcome of OVT to identify the optimal parameters for maximizing the efficacy of ZIKV–GSC virotherapy. Our numerical results showed that (1) the efficacy of the therapy was high as the SOX2 expression level was high, (2) two bifurcation values for both SOX2 expression level and viral bursting size were observed in GSC population dynamics, and (3) viral bursting size had symmetry related to the SOX2 expression level in terms of the stability of equilibrium points.

Materials and Methods

Model

2.1.

Experimental studies have shown that ZIKV-infected GSCs express the stem cell maker SOX2, with more than 90% of infected cells being SOX2^+^ [28]. The level of SOX2 expression correlates with susceptibility to ZIKV infection, highlighting SOX2 as a key determinant in the interaction between ZIKV and GSCs [27,29]. In this model, we assume that (1) GSCs proliferate with the logistic growth rate $[eqn]$ , up to their carrying capacity [27,41], (2) oncolytic ZIKV is 100% GSC-specific, (3) one virus particle infects one GSC; if a virus enters a GSC, it is incapable of infecting additional GSCs and ceases to be part of the free virus population, and (4) SOX2 expression level is included in the rate of cells infected by ZIKV. The three-dimensional model of GSCs with ZIKV is given by

$$[eqn]$$

$[eqn]$ is the population of GSCs; $[eqn]$ is the population of GSCs infected by ZIKV and the subpopulation of $[eqn]$ ; $[eqn]$ is the free Zika viruses;The term $[eqn]$ describes the logistic growth rate of GSCs; $[eqn]$ is the normalized rate of SOX2 expression level in a GSC, which ranges between 0 and 1;The constant value $[eqn]$ represents the strength of infectivity of the Zika virus in the GSCs;The term $[eqn]$ describes the rate of infected cells by free virus, $[eqn]$ ; $[eqn]$ is the bursting size of free virus particles; $[eqn]$ represents the death rate of infected GSCs after the cell oncolysis; $[eqn]$ is the clearance rate of the virus.

For non-dimensionalization, we set $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ . Then,

$$[eqn]$$

The system of Equation (1) becomes

$$[eqn]$$

We have the following model by setting the parameters:

$[eqn]$ , $[eqn]$ and $[eqn]$ . For convenience, we write $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ .

Then the system of Equation (3) becomes

$$[eqn]$$

Note that all the parameters ( $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ ) in the system of Equation (4) are positive constants.

Theorem 1. If $[eqn]$ , $[eqn]$ , and $[eqn]$ , then $[eqn]$ , $[eqn]$ , and $[eqn]$ for $[eqn]$ . Furthermore, the component $[eqn]$ is bounded and belongs to the interval [0, 1] for all $[eqn]$ .

Proof. If the conclusion $[eqn]$ , $[eqn]$ , and $[eqn]$ for $[eqn]$ is not true, then there exists a time $[eqn]$ at which at least one component first becomes zero. We examine each possible case.

Case 1: If the one component is zero at $[eqn]$ .

If $[eqn]$ , then $[eqn]$ . From the first equation of the system of Equation (4), $[eqn]$ for all $[eqn]$ by the uniqueness of the solution. Then, the second equation of the system of Equation (4) becomes $[eqn]$ . Using the separation of variable method, $[eqn]$ . This implies $[eqn]$ , since $[eqn]$ . From the third equation of the system of Equation (4), we have $[eqn]$ . This is a non-homogeneous linear differential equation because of the term by. The homogeneous solution is $[eqn]$ and the particular solution is $[eqn]$ . The general solution is $[eqn]$ . Therefore, $[eqn]$ .If $[eqn]$ , then $[eqn]$ . The first equation of the system of equation (4) is $[eqn]$ . Using the separation of variable method over $[eqn]$ , $[eqn]$ . So, $[eqn]$ for all $[eqn]$ . From the third equation of the system of Equation (4), we have $[eqn]$ . Similarly, $[eqn]$ . Thus, $[eqn]$ . Therefore, $[eqn]$ when $[eqn]$ .If $[eqn]$ , then $[eqn]$ . So, $[eqn]$ for $[eqn]$ , since $[eqn]$ , $[eqn]$ for all $[eqn]$ .

Case 2: If the two components are zero simultaneously at $[eqn]$ , it is easy to show that the third component will be non-negative for all $[eqn]$ .

Case 3: If the three components are zero simultaneously at $[eqn]$ , then $[eqn]$ , $[eqn]$ , and $[eqn]$ for all $[eqn]$ by the uniqueness of the solution.

Therefore, $[eqn]$ , $[eqn]$ , and $[eqn]$ for all $[eqn]$ . □

For boundness, from the first equation of the system of Equation (4), $[eqn]$ . Let $[eqn]$ , with initial condition $[eqn]$ , then the solution of the differential equation can be obtained by the separation of variable method and is given by $[eqn]$ . Since $[eqn]$ and $[eqn]$ , we have $[eqn]$ .

Therefore, all the solutions of the system of Equation (4) are non-negative and $[eqn]$ is bounded in the following region:

$$[eqn]$$

$[eqn]$ is the positive invariant where every solution with initial condition in $[eqn]$ remains there for all $[eqn]$ .

Analysis and Stability of Equilibrium

2.2.

The equilibrium points of the system are obtained by setting the right-hand side of the system of Equation (4) to zero. Let $[eqn]$ and $[eqn]$ .

Then, the system can be written as the autonomous system $[eqn]$ . We assume that the solution set $[eqn]$ of the system of Equation (4) are in $[eqn]$ . Simply, the equilibrium points are the solution of $[eqn]$ or $[eqn]$ which is given by

$$[eqn]$$

If $[eqn]$ , then $[eqn]$ and $[eqn]$ from the second and the third equations in Equation (5). Therefore, we have an equilibrium point $[eqn]$ .If $[eqn]$ and $[eqn]$ , we get $[eqn]$ from the second equation, which results in $[eqn]$ from the first equation in Equation (5). Then, $[eqn]$ . Thus, we have an equilibrium point $[eqn]$ .If $[eqn]$ and $[eqn]$ , then, from the second and third equations in Equation (5), $[eqn]$ . Since $[eqn]$ , we have $[eqn]$ . From the first and second equation in Equation (5), $[eqn]$ . Then, $[eqn]$ . From the third equation in Equation (5), $[eqn]$ . Then, $[eqn]$ . Thus, we have an equilibrium point $[eqn]$ , where $[eqn]$ , $[eqn]$ , and $[eqn]$ .

Basic Reproduction Number

2.3.

We consider the basic reproduction number as the number of secondary cases of infection generated from a single Zika virus in a GSC population where all GSCs are susceptible to infection. We utilize the next-generation matrix technique to calculate the basic reproductive number for the system of Equation (4).

Let $[eqn]$ , then all three equations of the system of Equation (4) can be written as $[eqn]$ , where

$$[eqn]$$

The Jacobian matrix of $[eqn]$ and $[eqn]$ at $[eqn]$ is

$$[eqn]$$

The basic reproduction number is obtained from the eigenvalue of $[eqn]$ , so we have

$$[eqn]$$

Stability of Equilibrium Points

2.4.

We studied the local stability of the equilibrium points using the linear stability analysis by finding the eigenvalues of Jacobian matrix at each equilibrium point. The Jacobian matrix of the nonlinear system in Equation (5) is given by

$$[eqn]$$

Theorem 2. The equilibrium point $[eqn]$ is always unstable.

Proof. The Jacobian matrix at $[eqn]$ is

$$[eqn]$$

The eigenvalues are $[eqn]$ , $[eqn]$ and $[eqn]$ . All parameters are positive values. $[eqn]$ and $[eqn]$ are negative but $[eqn]$ is positive. Therefore, the equilibrium point $[eqn]$ is unstable. □

Theorem 3. The equilibrium point $[eqn]$ is locally asymptotically stable if $[eqn]$ . Otherwise, it is unstable.

Proof. The Jacobian matrix at $[eqn]$ is given by

$$[eqn]$$

Solving the characteristic equation, $[eqn]$ , we have $[eqn]$ since $[eqn]$ and $[eqn]$ and $[eqn]$ satisfy the equation $[eqn]$ , where $[eqn]$ and $[eqn]$ . By Routh–Hurwitz criteria [42], the eigenvalues $[eqn]$ and $[eqn]$ are negative if $[eqn]$ and $[eqn]$ . Clearly, $[eqn]$ since all parameters are positive and for $[eqn]$ to be positive, $[eqn]$ , which implies $[eqn]$ . Therefore, if $[eqn]$ , $[eqn]$ is asymptotically stable. Otherwise, $[eqn]$ becomes unstable. Furthermore, $[eqn]$ , which leads $[eqn]$ . □

The equilibrium point $[eqn]$ represents a scenario where the population of GSC is at the carrying capacity (or maximum value of 1), indicating no infection has taken place. This means that the virotherapy has failed completely, as all GSCs remain uninfected. The theorem states that $[eqn]$ is locally asymptotically stable if $[eqn]$ . This implies that when the basic reproduction number $[eqn]$ is less than 1, each infected cell generates, on average, less than one new infected cell, leading to the decline and eventual extinction of the infection. As a result, the therapy is infective and the GSCs remain intact. Conversely, if $[eqn]$ is greater than or equal to 1, the infection can spread, making $[eqn]$ unstable.

Theorem 4. The equilibrium point $[eqn]$ is locally asymptotically stable if $[eqn]$ and $[eqn]$ .

Proof. The Jacobian matrix at $[eqn]$ is

$$[eqn]$$

Since $[eqn]$ from the first equation in the system of Equation (5),

$$[eqn]$$

The characteristic equation is $[eqn]$ , □

where $[eqn]$ , $[eqn]$ , and $[eqn]$ . Therefore, we have $[eqn]$ , $[eqn]$ , and $[eqn]$ . It is clear that $[eqn]$ and $[eqn]$ since the equilibrium point $[eqn]$ exists if $[eqn]$ and $[eqn]$ . It is enough to show that $[eqn]$ to prove that the eigenvalues of $[eqn]$ are negative by Routh–Hurwitz criteria.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Now, the term $[eqn]$ is always positive, so we can ignore it when considering the sign of $[eqn]$ . Therefore, $[eqn]$ holds if and only if $[eqn]$ . This inequality is equivalent to $[eqn]$ . Therefore, we have shown that $[eqn]$ if and only if $[eqn]$ . Note that $[eqn]$ .

Therefore, from Theorems 3 and 4, the stability of $[eqn]$ and $[eqn]$ changes about $[eqn]$ , which means a transcriptical bifurcation happens at $[eqn]$ .

The equilibrium point $[eqn]$ corresponds to a situation where there is a balance between the populations of GSCs, infected GSCs, and the ZIKV. This balance indicates that the infection is present and affects GSCs. The theorem states that $[eqn]$ is locally asymptotically stable if $[eqn]$ . This implies that when the basic reproduction number $[eqn]$ is greater than 1, each infected GSC generates more than one new infected GSC, leading to the sustained spread of the infection. In this scenario, the ZIKV effectively infects the GSCs and reduces their population indicating successful virotherapy.

Estimation of Parameters

2.5.

Parameters were estimated from the experimental data (Table 1) [28]. In the experiment, there were 1 million cancer stem cells on day 0, and the number of cancer stem cells increased to 9 million by day 7. The difference in cell count from day 0 to day 7 was 8 million cells, and the population followed exponential growth. $[eqn]$ . In the experiment [28], they infected GSCs with ZIKV at a multiplicity of infection (MOI) of 5. More than 60% of GSCs were infected 48 h after infection, and 90% of the infected cells expressed SOX2. We estimated the infection rate as follows: the initial number of GSCs was 1 million, and the MOI was 5. Thus, the number of infectious viral particles was 1 million $[eqn]$ million. Therefore, the number of infected cells at 48 h was 60% $[eqn]$ million $[eqn]$ million. Given that 90% of the infected cells expressed SOX2, the number of SOX2-expressing infected cells was 90% $[eqn]$ million = 2.7 million. The infection rate can be estimated as the fraction of infected cells produced per unit time per initial number of cells. Assuming exponential growth of infection over time, we estimated the average rate of increase in the number of infected cells. Since we observed 3 million infected cells at 48 h, we can estimate the average rate of increase over the first 48 h. Therefore, the infection rate would be 2.7 million/(1 million $[eqn]$ days) = 0.054 per day. The ZIKVs were cleaned up after day 8, so the death rate of ZIKV was 1/8 = 0.125 per day.

Sensitivity Analysis

2.6.

Sensitivity analysis is applied to study the effect of parameters on the proposed mathematical model. In particular, it is necessary to identify the most sensitive parameters that cause a disturbance in the model dynamics with a small change in their numeric values. To check the sensitivity of $[eqn]$ , sensitivity, we calculate its derivatives as follows:

$$[eqn]$$

Since all the parameters are positive, $[eqn]$ , $[eqn]$ , and $[eqn]$ . It concludes that the basic reproduction number $[eqn]$ increases as $[eqn]$ , $[eqn]$ and $[eqn]$ increase. The normalized sensitivity indices corresponding to these parameters are estimated as follows:

$$[eqn]$$

Here, the sensitivity index can be constant depending on some parameters or can be free of any independent parameters. The partial rank correlation coefficient (PRCC) results for significance of parameters involved in $[eqn]$ is shown in Figure 1. The positive PRCC values for $[eqn]$ , $[eqn]$ , and $[eqn]$ indicate that these parameters are directly correlated with $[eqn]$ , meaning that increases in these parameters lead to higher $[eqn]$ values, thereby enhancing the infection spread. On the other hand, the negative PRCC value for $[eqn]$ demonstrates its inverse relationship with $[eqn]$ ; as $[eqn]$ increases, $[eqn]$ decreases, suggesting that enhancing virus clearance is crucial for controlling the infection. This sensitivity analysis emphasizes the importance of accurately estimating these parameters, as small changes in their values can significantly alter the infection dynamics, guiding effective intervention strategies.

Results

The nondimensionalized model (Equation (4)) was employed to present the numerical results. For these calculations, we utilized the Runge–Kutta 2nd order method with a time step of $[eqn]$ in MATLAB R2023a (The MathWorks, Natick, MA, USA). To assess the accuracy of the numerical scheme, we also tested smaller values of $[eqn]$ and compared the results with those obtained using the Runge–Kutta 4th order method.

Existence and Stability of the Equilibrium Points

3.1.

For numerical simulations, we initially set parameters as $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ . However, we may adjust certain parameter values to ensure that the existence and stability conditions for each equilibrium point are met. Figure 2 illustrates the population solutions over relative time (A and C) and the trajectories of solutions in the phase space (B and D). We set the SOX2 expression level constant to the value $[eqn]$ , which satisfies the condition $[eqn]$ , ensuring that $[eqn]$ becomes asymptotically stable (Figure 2A,B). The equilibrium point $[eqn]$ becomes asymptotically stable when $[eqn]$ , while $[eqn]$ becomes asymptotically unstable and it satisfies the condition $[eqn]$ . The solution $[eqn]$ converges to $[eqn]$ . This result provides a good agreement between our analytic and numerical results.

SOX2 Expression Level Changes the Structure of GSC Dynamics

3.2.

The bifurcation diagram of GSCs with respect to the normalized SOX2 expression levels is shown in Figure 3. In this simulation, we set the virus bursting size $[eqn]$ . Our numerical result shows two bifurcation values: (1) A transcritical bifurcation occurs at $[eqn]$ (green colored vertical line in Figure 3) which results in a qualitative change in stability between equilibrium points $[eqn]$ and $[eqn]$ . When $[eqn]$ , $[eqn]$ is stable, while $[eqn]$ is unstable. However, the stability of two equilibrium points changes for $[eqn]$ ; $[eqn]$ becomes stable and $[eqn]$ becomes unstable. (2) A Hopf bifurcation is observed at $[eqn]$ (red colored vertical line in Figure 3), which leads to oscillatory behavior in populations. When $[eqn]$ , $[eqn]$ is a stable equilibrium point. However, when $[eqn]$ , the system undergoes a qualitative change in dynamics such as a periodic cycle around $[eqn]$ . This result explains the sensitivity of the system to SOX2 expression level, with $[eqn]$ which represents a critical value where periodic patterns become prominent in the GSC–ZIKV dynamical system.

From a biological point of view, the dynamic interplay among GSCs, infected GSCs, and ZIKV has direct implications for therapeutic efficacy. For example, equilibrium point $[eqn]$ represents a free virus equilibrium point where therapy fails (GSC population approaches to the carrying capacity). The other equilibrium point $[eqn]$ represents partial success which indicates a reduction in the GSC population. This bifurcation analysis provides a good understanding of the complex population dynamics within the GSC–ZIKV system. These results suggest a critical threshold of SOX2 expression level is an important factor to be considered for therapeutic strategies to eradicate the GSCs in the tumor microenvironment.

Interplay between SOX2 Expression Level and Bursting Size Affects the Dynamic of OVT

3.3.

For different values of the bursting size (b), we investigated the effect of SOX2 expression level (s) on GSC dynamics with different values of the bursting size value $[eqn]$ shown in Figure 4 (A–D were when $[eqn]$ , 5, 15, and 20, respectively). For $[eqn]$ , the GSC population converged to the carrying capacity, which resulted in the therapy failure. For $[eqn]$ , a transcritical bifurcation occurred at $[eqn]$ where $[eqn]$ became stable while $[eqn]$ became unstable. The GSC population reached its minimum of 0.78 at $[eqn]$ . Increasing the bursting size to 15 resulted in shifting a transcritical bifurcation threshold to $[eqn]$ . The minimum GSC population, 0.23, was observed at $[eqn]$ . Finally, for $[eqn]$ , there were two bifurcations threshold values: a transcritical bifurcation at $[eqn]$ and a Hopf bifurcation at $[eqn]$ . At $[eqn]$ , the minimum population was 0.167. Moreover, when $[eqn]$ , all populations showed oscillations. These results indicate a shift in the transcritical bifurcation point $[eqn]$ to a lower value of the bursting size $[eqn]$ and the induction of an oscillatory pattern with higher bursting size.

Both SOX2 expression level (s) and virus bursting size (b) play an important role in OVT since (1) SOX2 expression levels activate the membrane receptor αvβ5, which results in an enhancing of the infection rate of ZIKV into GSCs and (2) the virus bursting size is directly proportional to the basic reproduction number $[eqn]$ , which can change the stability of equilibrium points and affect the efficacy of the therapy. Understanding the underlying mechanisms of changes of stability of equilibrium points or of structure of dynamics (or bifurcation) helps to show the importance of the optimization of parameters in the system and in treatment efficacy.

Stability Regions for Equilibrium Points in Two-Dimensional Parameters: SOX2 Expression Level and Bursting Size

3.4.

In our mathematical exploration, we systematically varied the bursting size (b) and SOX2 expression level (s), conducting simulations to calculate eigenvalues and visualize stability regions for different equilibrium points on the s-b coordinate. The two-dimensional bifurcation diagram is shown in Figure 5. The color-coded representation denotes blue for the stability region of equilibrium point $[eqn]$ , green for $[eqn]$ stability, and yellow for regions where oscillations occur. The boundary between the green and yellow regions highlights the minimum GSC population. These stability regions elucidate critical points in the interplay between SOX2 expression, bursting size, and GSC dynamics. Transitions among the blue, green, and yellow regions signify qualitative shifts, reflecting the impact of SOX2 expression levels on GSCs. Understanding these stability regions is crucial from a therapeutic perspective. The blue region suggests potential resistance scenarios, while the green and yellow regions offer opportunities for targeted interventions. The boundary representing the minimum population holds significance, indicating conditions conducive to maximizing therapeutic efficacy. The symmetrical patterns observed in the bifurcation diagram provide a framework for predicting the behavior of the system. The interplay between SOX2 expression level and the virus bursting-size parameters reveals a symmetric pattern in the stability region of equilibrium points.

Discussion

In this study, we established a simple ODE mathematical model to understand the interactions among GSCs, infected GSCs, and ZIKV with a primary focus on the role of SOX2 expression level and the bursting size of ZIKV in dynamic populations. We provided analytic work on the existence and boundness of equilibrium points or solutions. The stability of equilibrium points was performed by local stability analysis. As a result, we showed that the stability of equilibrium points is dependent on the basic reproduction number $[eqn]$ , such as if (1) $[eqn]$ is stable with $[eqn]$ , and (2) $[eqn]$ becomes stable with $[eqn]$ . Sensitivity analysis was evaluated, and the Runge–Kutta 2nd order method was used for numerical simulations.

In our experimental results, we identified two threshold values of SOX2 expression level where SOX2 expression can change from (1) free ZIKV equilibrium point to equilibrium point with ZIKV (transcritical bifurcation) and (2) steady state solutions to oscillations (Hopf bifurcation). We explored how bursting size and SOX2 expression level relate to the efficacy of virotherapy by conducting simulations to analyze stability regions for various equilibrium points on a coordinate system. The effect of SOX2 expression level and bursting size of ZIKV on the stability of equilibrium points is shown in the 2D bifurcation diagram (Figure 5). These two parameters play a crucial role in the dynamics of OVT for therapeutic strategies. The boundary separating green and yellow zones represents the minimum GSC population, offering crucial insights for virotherapy. Transitions between these regions reflect qualitative shifts influenced by SOX2 expression and busting size. Understanding these stability regions is essential for developing effective virotherapy strategies and optimizing therapeutic outcomes.

Mathematical modeling serves as a powerful tool in the development and optimization of therapeutic strategies against many human diseases, including brain tumors [43–45]. Moreover, it has the potential to provide new insights, hypotheses, and experimental directions, ultimately leading to the development of personalized cancer therapies [43–45]. Our study aimed to develop a mathematical model to precisely predict the optimal levels of SOX2 and bursting size for achieving successful virotherapy outcomes against GSCs. The high recurrence and therapy resistance of human glioblastoma may be attributed to the presence of GSCs in GBM, which often accompany high expression levels of the transcription factor SOX2. The Zika viruses are known to target GSCs via SOX2-integrin-mediated infection [27,29]. Therefore, SOX2 can be used as a prediction marker for successful virotherapy as the expression level of SOX2 in GSCs can predict their susceptibility to ZIKV infection and subsequent virotherapy response.

Although our study focused on the aspect of the role of SOX2 expression levels and ZIKV infection, the limitations may include the exclusion of immune responses in our modeling approaches. During oncolytic virotherapy, immune responses, including the activation of natural killer (NK) cells, may influence the anti-tumor efficacy of virotherapy [37]. Our previous study established a mathematical model to assess the fine balance between viral replication rates and NK cell activities during virotherapy. Additionally, recent research has shown that ZIKV E protein co-localizes with SOX2, generating a long-term memory for antitumor immune response [46]. Conducted on mice, this study showed that mice treated with a combination of ZIKV lived significantly longer than other groups, surviving for more than 120 days after treatment [46]. Moreover, ZIKV treatment not only elevated cytotoxic T cell infiltration but also enhanced the antitumor immune response within GSCs by recruiting and activating T cells [46]. In a recent study by Garcia et al., the attenuated Zika viral strain ZOL-1 was utilized in a mouse model of GBM [47]. Their findings underscore the importance of molecular characterization of the target cells for successful virotherapy, as they found two different responsive cancer groups, showing that some patient-derived cancer cells responded better than others [47]. These studies emphasize the necessity for further investigation into the correlation between immune responses and viral infection to achieve maximum efficacy in ZIKV-mediated virotherapy.

Conclusions

Our mathematical model determined the optimal SOX2 levels for efficient and enhanced virotherapy by identifying critical threshold levels of both SOX2 and viral replication. Given that SOX2 levels in human glioblastoma stem cells are critical for successful ZIKV infection and replication, our model enables precise prediction of the success rates of ZIKV-targeted GSC lysis. This study highlights the potential for developing personalized virotherapies based on a patient’s SOX2 expression levels in their GSCs. Since immune responses may significantly impact the efficacy and outcomes of oncolytic virotherapy, future studies should focus on elucidating the optical conditions for recruiting immune cells, such as cytotoxic T cells, to the infected GSCs and glioblastoma cells to maximize the effects of virotherapy against GBM.