Sensitivity analysis unveils the interplay of drug-sensitive and drug-resistant Glioma cells: Implications of chemotherapy and anti-angiogenic therapy

1.   Introduction

Gliomas represents the most prevalent category of primary brain tumors, encompassing their exceptionally aggressive variant known as glioblastoma multiforme (GBM). This subtype constitutes approximately 15% of all brain tumors ^[1]. Malignant gliomas are aggressive brain tumors known for their rapid development of blood vessels (angiogenesis), which is essential for their growth in the brain. These tumors exhibit a high level of proliferation of endothelial cells, a key feature in their classification according to the World Health Organization (WHO) classification system. This process of angiogenesis, which involves complex interactions between tumor cells and blood vessel cells, plays a critical role in the behavior of these tumors and the prognosis of patients ^[2].

Chemotherapy remains a potential approach for treating cancer despite these advances. Failures in chemotherapy are associated with drug resistance. Drug resistance is now a major problem in the field of cancer. The long-term efficacy of drugs aimed at cancer patients is frequently inevitably constrained by drug resistance. Thousands of efforts have been put toward reducing drug resistance and increasing patient survival ^[3].

Chemotherapy combinations often include antiangiogenic treatment. Antiangiogenic treatment is a technique used to treat cancer to obstruct the ability of blood vessels to carry nutrients and oxygen to tumor cells while also halting the growth of new blood vessels. Vascular endothelial growth factor (VEGF), which is considered a primary promoter of angiogenesis, is the target of the majority of licensed antiangiogenic drugs used to treat cancer ^[4,5]. Because VEGF also exhibits immunosuppressive properties, which emphasizes a potential target for antiangiogenic therapy, these medications can improve immunotherapy in addition to reducing angiogenesis. Given the high rate of endothelial proliferation, increased vascular permeability, and the production of proangiogenic growth factors ^[6], targeting blood vessels in brain tumors has become a very attractive method.

A relatively recent concept in cancer research is the idea of cancer dormancy, which is an additional characteristic hallmark of cancer. Angiogenic and immunogenic dormancy processes are well known; there is also a form of cellular dormancy that occurs within the tumor at the individual cell level. Before angiogenesis undergoes changes, increased cell growth leads to a decrease in oxygen and nutrient levels in areas far from blood vessels. This, in turn, causes cell death and establishes a balance between cancer cell growth and cell death ^[7,8]. Dormant tumor cells are quite common in the general population ^[9], and those that persist after primary tumor treatment or removal often exhibit resistance to chemotherapy ^[10,11]. Predictions in the case of mathematical modeling are helpful to understand the extent of the disease population, the effects of widespread disease, and how long the disease will last. The author also performed a sensitivity analysis of the model to identify which parameters had an impact and to interpret their biological significance. By analyzing the effect and extent of the spread of the disease through a mathematical model, accurate predictions can be made regarding how the infection will spread so that preventive and treatment measures can be taken against the spread of the disease in a population.

2.   Model description and analysis

This model is based on previous research by ^[12], incorporating two types of disease cells: Those sensitive to drugs and those resistant to drugs. It is important to note that in this model, we considered antiangiogenic therapy as a form of continuous treatment. Within this model, we also account for dormancy occurring within cells, called angiogenic dormancy. This factor influences the balance between cell proliferation and cell death.

These types include cells affected by the disease, namely the concentration of glioma-sensitive cells ($g_2$) and glioma-resistant cells ($g_3$), with healthy cells among them, such as the concentration of glial cells ($g_1$), endothelial cells ($g_4$), and neurons ($g_5$), as well as the concentrations of chemotherapy ($q$) and antiangiogenic agents ($y$). The parameters of the non-dimensionalized model used in the model are shown in Table 1:

As per the reference ^[12], the non-dimensionalized model can be expressed as follows:

with

and initial values $g_i\geq0, i = 1, ..., 5$ $q\geq0, y\geq0$ for $t = 0$ and $F(x)$ is a function defined as

2.1.   Basic reproduction number

In this subsection, we calculate the basic reproduction numbers for a model involving a disease with drug-sensitive and drug-resistant strains. The basic reproduction number (${{R}_{0}}$) measures the rate of spread of a tumor. If ${{R}_{0}}$ is greater than one, it suggests that the number of affected cells, which encompasses both drug-sensitive and drug-resistant cases, will increase, implying the persistence of the affected cells. On the other hand, if is less than one, it indicates that, on average, each tumor cell produces fewer than one new cell and, therefore, therapy (administration of drugs) has the potential to eradicate the tumor. In this model, at each time step, a tumor cell gives rise to offspring or dies, serving as a parameter that determines whether the tumor will continue growing or will be suppressed and eventually eliminated by therapy ^[13,14]. In our analysis, we employ the next-generation matrix technique to estimate the basic reproduction numbers for our system, which encompasses both drug-sensitive and drug-resistant forms of the disease ^[15].

The simplified model considers two disease states: Drug-sensitive (${{g}_{2}}$) and drug-resistant (${{g}_{3}}$) states, along with five non-disease states: Glial cell (${{g}_{1}}$), endothelial cells (${{g}_{4}}$), and neurons (${{g}_{5}}$), as well as chemotherapy agent ($q$) and antiangiogenic agent ($y$). According to ^[12], to determine the glioma-free equilibrium point, we evaluate the system of Eqs (2.1)–(2.7) when $g_2 = g_3 = 0$. The first glioma-free equilibrium point is $E_0 = (0, 0, 0, 0, 0, q, y)$, which represents the non-cell state. Therefore, from Eqs (2.6) and (2.7), in the absence of glioma, we have $q = \frac{\phi}{\psi}$ and $y = \frac{\delta}{\gamma}$.

When we linearize the system around this first glioma-free equilibrium, we discover that the Eqs (2.2) and (2.3) govern the dynamics of and create a closed system, resulting in a linearized sub-model for the disease dynamics involving both drug-sensitive and drug-resistant strains. Let $x$ represent the disease compartments and $y$ represent the non-disease compartments, so it can be written as follows,

and

so that $\dot{{{\bf x}}} = \mathcal{F}({{\bf x}}, {{\bf y}})-\nu({{\bf x}}, {{\bf y}})$ where

The matrix $\mathcal{F}$ contains the transmission component of ${{\bf x}}$ (i.e., the arrival of susceptible individuals into the disease compartments $g_2$ and $g_3$) and the matrix $\nu$ contains transitions between, and out of the disease states (i.e., mutation, dormancy and death), then $\dot{{{\bf y}}} = {{\bf s(x, y)}}$ is as follows:

Before calculating the basic reproduction number using the next-generation matrix, several assumptions need to be met, as follows:

1) Based on $\mathcal{F}({{\bf x}}, {{\bf y}})$ and $v({{\bf x}}, {{\bf y}})$ in Eq (2.9), it is obtained that $\mathcal{F}({\bf 0}, {{\bf y}}) = 0$ and $v({\bf 0}, {{\bf y}}) = 0$ for every ${{\bf y}} \geq {\bf 0}$.

2) Based on Eq (2.10), the disease-free system ${\bf \dot{y}} = {\bf s}({\bf 0}, {{\bf y}})$ has a stable asymptotic equilibrium point. The disease-free population in the system (2.1)–(2.7) is $g_1, g_4, g_5, q$, and $y$, so from Eq (2.10), we have:

The first equilibrium point for glioma from system ${\bf \dot{y}} = {\bf s}({\bf 0}, {{\bf y}})$ is $(0, 0, 0, \frac{\phi}{\psi}, \frac{\delta}{\gamma})$. The Jacobian matrix of the system ${\bf s}({\bf 0}, {{\bf y}})$ at the point $(0, 0, 0, \frac{\phi}{\psi}, \frac{\delta}{\gamma})$ is as follows:

Next, we calculate the eigenvalues of the matrix $D f({\bf s}({\bf 0}, {{\bf y}}))$ as follows:

The eigenvalues obtained are $\lambda_1 = p_1 - \left(d_{12}\frac{\delta}{\gamma} + d_{10}\right)\frac{\phi}{\psi a_1}$, $\lambda_2 = p_4-\frac{d_4 \delta}{\gamma a_4}$, $\lambda_3 = 0$, $\lambda_4 = -\psi$, and $\lambda_5 = -\gamma$. Since $\psi > 0$ and $\gamma > 0$, we have $\lambda_4 < 0$ and $\lambda_5 < 0$. Furthermore, if $\phi > \frac{p_1\psi a_1 \gamma }{(d_{10} \gamma+d_{12} \delta)}$, then $p_1 - \left(d_{12}\frac{\delta}{\gamma} + d_{10}\right)\frac{\phi}{\psi a_1} < 0$, and we obtain $\lambda_1 < 0$. Similarly, if $\delta > \frac{p_4\gamma a_4 \delta}{d_4}$, then $p_4-\frac{d_4}{\gamma a_4} < 0$, and we obtain $\lambda_2 < 0$. Therefore, as for all $\lambda_i < 0$ for $i = 1, 2, 3, 4, 5$, the equilibrium point of the disease-free system is asymptotically stable.

3) The values of $\mathcal{F}({{\bf x}}, {{\bf y}}) \geq {\bf 0}$ for every ${{\bf x}}, {{\bf y}} \geq {\bf 0}$.

4) If ${{\bf x}} = {\bf 0}$, then $v({{\bf x}}, {{\bf y}}) \leq {\bf 0}$.

5) Based on the matrix $v({{\bf x}}, {{\bf y}})$, it is obtained that

so that $\sum_{i = 1}^2 v_i({{\bf x}}, {{\bf y}}) \geq 0$ for every ${{\bf x}}, {{\bf y}} \geq {\bf 0}$.

The system (2.1)–(2.7) satisfies the five assumptions of the next-generation matrix. Therefore, to calculate the basic reproduction number, the next-generation matrix method can be used. Next, we will determine the matrix $F$, which is the Jacobian matrix of the matrix $\mathcal{F}$ at the disease-free equilibrium point, and then the matrix $V$ will be the Jacobian matrix of the matrix $v$ at the first equilibrium point of glioma. Here is the matrix $F$ and $V$ at the first glioma-free equilibrium point:

The next-generation matrix, $M$, is then given by ^[15]

The eigenvalues of $M$ obtained are $\lambda_1 = \frac{{p_{2} \gamma \psi a_{2}}}{{(a_{2} (u + \rho) \psi + d_{20} \phi) \gamma + d_{22} \delta \phi}}$ and $\lambda_2 = \frac{{p_3}}{{\rho}}$. The principal eigenvalues of matrix $M$ serve as the fundamental reproduction rates for both the drug-susceptible and drug-resistant glioma. They represent the average number of new infections generated by a single infected individual from each strain. The lower triangular structure of $M$ allows an immediate extraction of the fundamental reproduction rates for the drug-susceptible and drug-resistant glioma, respectively, as follows:

and

so that

At $R_{01}$, we choose either $R_{0 A}$ or $R_{0 B}$ by taking the maximum between them.

Analogously to the steps in the first reproduction number, we obtain the second reproduction number associated with the disease-free state $E_1 = (g_1^b, 0, 0, g_4^b, 0, q^b, y^b)$. Therefore, we have the second reproduction number

so that

At $R_{02}$, we choose either $R_{0 C}$ or $R_{0 D}$ by taking the maximum between them. Interestingly we find that the basic reproduction numbers $g_2$ and $g_3$ are both independent of the amplification rate $\rho$ ^[16].

2.2.   System properties

2.2.1.   Existence equilibria

Equations (2.1)–(2.7) indicate the existence of a glioma-free equilibrium, $E_0$, as shown in ^[12].

It is also clear that there exists a second glioma-free equilibrium, denoted as ${{E}_{1}}$, which is always present and defined as $E_1 = (g_1^b, 0, 0, g_4^b, 0, q^b, y^b)$.

For the existence of the equilibriua $E_1$ note that if $\delta\left(d_4\right) < a_4 \gamma p_4$, then $-\gamma a_3+\gamma+c_4 < \left(\gamma a_4-\gamma-c_4\right)^2 - 4\left(\gamma+c_4\right)\left[\delta\left(d_4\right) / p_4-a_4 \gamma \right]$. Therefore, if $\phi/\psi < p_1 a_1/(d_{10}+d_{11} g_4^b+d_{12} y^b)$, $g_1^b > 0$ always exists, as well as for $q^b$.

Furthermore, from Eqs (2.1)–(2.7) we can also derive the existence of a mono-existent endemic equilibrium, denoted as ${{E}_{2}}$, where the drug-resistant strain persists while the drug-susceptible strain decreases ${{E}_{2}} = (0, 0, g_{3}^{r}, g_{42}^{r}, 0, {{q}^{r}}, {{y}^{r}})$, where

and $g_{42}^r$ is a real positive solution of the following equation: ${{l}_{1}}g_{4}^{3}+{{l}_{2}}g_{4}^{2}+{{l}_{3}}{{g}_{4}}+{{l}_{4}} = 0$, where $l_i$, for $i = 1, 2, 3, 4$, are defined as follows:

When we examine Eq (2.11), we can observe that the mono-existence of the endemic equilibrium exists if and only if $R_{0B} \geq 1$, $3 l_3 < l_2^2$, $\frac{q^2}{4} > -\frac{p^3}{27}$, $2 l_2^3+27 l_4 < 9 l_2 l_3$ and $l_2 < 0$.

2.3.   Stability analysis

In this subsection, the stability of each equilibrium point of the system (2.1)–(2.7) will be analyzed through linearization. The following results are established:

Lemma 1. If ${{R}_{01}} = \max [{{R}_{0A}}, {{R}_{0B}}] < 1$, $\phi > \frac{p_1 \psi a_1 \gamma}{(d_{10} + d_{12} \delta)}$, $\delta > \frac{p_4 \gamma a_4}{d_4}$ the non-cell equilibrium $E_0 = (0, 0, 0, 0, 0, \frac{\phi}{\varphi}, \frac{\delta}{\gamma})$ is locally asymtotically stable: If, however, ${{R}_{01}} = \max [{{R}_{0A}}, {{R}_{0B}}] > 1,$ at least one of the eigenvalues has a positive real part, rendering ${{E}_{0}}$ unstable.

Proof. We consider the Jacobian of the system (2.1)–(2.7) at the first glioma-free equilibrium point, ${{E}_{0}}$, which reduces to which is given by $D{\bf f}E_{0}$. By symbolizing each component of the Jacobian matrix with $m_{ij}$, where $i$ represents the row index and $j$ represents the column index, we obtain

The structure of $D{\bf f}{E_0}$ allows us to immediately read off the seven eigenvalues, $\lambda_i$, as

It can be readily confirmed that for ${{R}_{0A}} < 1$, ${{R}_{0B}} < 1$, $\phi > \frac{p_1 \psi a_1 \gamma}{(d_{10} + d_{12} \delta)}$, and $\delta > \frac{p_4 \gamma a_4}{d_4}$ all eigenvalues have negative real parts. Consequently, the equilibrium of the non-cell ${{E}_{0}}$ in the system (2.1)–(2.7) is locally asymptotically stable under these conditions. However, if ${{R}_{0A}} > 1$ or ${{R}_{0B}} > 1$, at least one of the seven eigenvalues has a positive real part, making $E_0$ unstable.

Lemma 2. If ${{R}_{02}} = \max [{{R}_{0C}}, {{R}_{0D}}] < 1$, $i_1 > 0$, ${{i}_{1}}{{i}_{2}} > {{i}_{3}}$, $i_1 i_2 i_3 > i_3^2+i_1^2 i_4$ the first glioma-free equilibrium is locally asymptotically stable: If, however, ${{R}_{02}} = \max [{{R}_{0C}}, {{R}_{0D}}] > 1$, at least one of the eigenvalues has a positive real part, making ${{E}_{1}}$ unstable.

Proof. We consider the Jacobian of the system (2.1)–(2.7) at the glioma-free equilibrium point, ${{E}_{1}}$, which reduces to which is given by $D{\bf f}{E_1}$. By symbolizing each component of the matrix with $e_{ij}$, where $i$ represents the row index and $j$ represents the column index, we obtain

The structure of $D{\bf f}{E_1}$ allows us to immediately read off the first to third eigenvalues,

then we find $\lambda_4, \lambda_5, \lambda_6, \lambda_7$ from the roots of the following equation

with

For local stability, we must ensure that the Routh-Hurwitz criteria is satisfied; it will be negative if $i_1 > 0$, ${{i}_{1}}{{i}_{2}} > {{i}_{3}}$, $i_1 i_2 i_3 > i_3^2+i_1^2 i_4$. It can be easily confirmed that for ${{R}_{0C}} < 1$, ${{R}_{0D}} < 1$, all eigenvalues have negative real parts. Consequently, the second glioma-free equilibrium ${{E}_{1}}$ in the system (2.1)–(2.7) is locally asymptotically stable under these conditions. However, if ${{R}_{0C}} > 1$ or ${{R}_{0D}} > 1$, at least one of the seven eigenvalues has a positive real part, rendering ${{E}_{1}}$ unstable.

Lemma 3. Given the values of $s_2$, $B$, and $P$ in Eqs (2.15), (2.19), and (2.20), if $R_{0B} > max[1, R_{0A}]$, $p_1 < \beta_1g_3^r + \frac{(d_{11}g_4^r + d_{12}y^r + d_{10})q^r}{a_1}$, $\frac{P}{3 \sqrt[6]{2}} < \frac{\sqrt[2]{3} B}{3 P}+\frac{s_2}{3}$ and $\frac{P}{6 \sqrt[3]{2}}+\frac{s_2}{3} > \frac{\sqrt[3]{2} B}{6 P}$, then the equilibrium point $E_2 = (0, 0, g_3^r, g_4^r, 0, q^r, y^r)$ is asymptotically stable.

Proof. We consider the Jacobian of the system (2.1)–(2.7) at the glioma-free equilibrium point, ${{E}_{2}}$, which reduces to which is given by $D{\bf f}{E_2}$. By symbolizing each component of the Jacobian matrix with $o_{ij}$, where $i$ represents the row index and $j$ represents the column index, we obtain

The structure of $D{\bf f}{E_2}$ allows us to immediately read off the first to fourth eigenvalues,

then we find $\lambda_5, \lambda_6, \lambda_7, \lambda_8$ from the roots of the following equation

with

Next, the roots of the characteristic equation (2.14) are obtained by following the steps of Cardano's formula, provided by ^[17] as follows:

where

Next, we will analyze the real part of these eigenvalues and the conditions under which the real part of the eigenvalues is negative. Let

The condition for $P$ to be real is that $A^2 \geq B$ then we obtain:

The condition for $\lambda_5 < 0$ is $\frac{P}{3 \sqrt[6]{2}} < \frac{\sqrt[2]{3} B}{3 P}+\frac{s_2}{3}$, and the condition for $Re_{\left(\lambda_{6, 7}\right)} < 0$ is $\frac{P}{6 \sqrt[3]{2}}+\frac{s_2}{3} > \frac{\sqrt[3]{2} B}{6 P}$.

Since the real part of all eigenvalues $\lambda_i$ is negative for each $i = 1, 2, 3, 4, 5, 6, 7$, the equilibrium point $E_2 = (0, 0, g_3^r, g_4^r, 0, q^r, y^r)$ is asymptotically stable.

3.   Analysis of parameter sensitivity

Sensitivity analysis is performed to guide the parameters that contribute the most to cancer treatment efficacy. In this study, a normalization index method ^[18] will be employed to indicate which treatment parameters contribute the most to cancer eradication. The analysis is focused on parameters related to the basic reproduction number. A sensitivity analysis of this model is conducted to determine the impact of changes in parameter values on the values of the basic reproduction number. We focus on sensitivity analysis regarding the second reproduction number $R_{02}$, specifically the glioma-free state. If the data were available, we could use ordinary least squares and its extension^[30] to estimate the parameters. However, if the data were affected by uncertainties, the parameters could still be estimated using a fuzzy approach; see, e.g., ^[31,32,33] for a more detailed method. In this paper, however, we will only use predefined parameter values, as given in Table 2. The results of the sensitivity analysis for a parameter are as follows.

From Table 3, it is evident that the parameter with a positive sensitivity index is $\psi$. This indicates that if $\psi$ is increased while keeping the other parameters constant, it will increase the value of $R_{0C}$ and, consequently, increase the endemicity of tumor cells from glioma. On the other hand, the parameters $\phi$ and $\rho$ have negative values of the sensitivity index, meaning that if $\phi$ or $\rho$ increases while keeping the other parameters constant, it will decrease the value of $R_{0C}$ and consequently reduce the endemicity of glioma tumor cells.

From Table 4, it is evident that the sensitivity indices with positive values are associated with the parameters $\phi$. This indicates that if one of the parameters of $\phi$ is increased while keeping the other constant, it will increase the value of $R_{0D}$ and, consequently, increase the endemicity of tumor cells from gliomas. On the other hand, the parameters $\rho$ and $\psi$ have negative sensitivity indices. This means that if one of the parameters of is increased while keeping the others constant, it will decrease the value of $R_{0D}$ and, consequently, reduce the endemicity of glioma tumor cells.

In Figure 1, the sensitivity analysis indicates that the sensitivity index of the parameter $\rho$ to $R_{0C}$ is -0.2954. This means that an increase in 10% of the parameter $\rho$ while keeping other parameters constant will result in a decrease of 2.954% in $R_{0C}$. On the contrary, a 10% decrease in $\rho$ will lead to a 2.954% increase in $R_{0C}$. The analysis of sensitivity indicates that the $\rho$ parameter has a sensitivity index of -0.84746 in relation to $R_{0D}$. This suggests that if the parameter $\rho$ is increased by 10%, the value of $R_{0D}$ will decrease by 8. 4746%. Conversely, if $\rho$ is decreased by 10%, the value of $R_{0D}$ will increase by 8. 4746%.

In Figure 2, the sensitivity analysis results reveal that the sensitivity index of the parameter $\psi$ to $R_{0C}$ is 0.61836. This means that if the parameter $\psi$ increases by 10%, the value of $R_{0C}$ also increases by 6.1836%. Conversely, if the parameter $\psi$ decreases by 10%, the value of $R_{0C}$ decreases by 6.1836%. According to the sensitivity analysis findings, the sensitivity index of the parameter $\psi$ to $R_{0D}$ is approximately $-1.478 \times 10^{-7}$. In simpler terms, if the parameter $\psi$ is increased by 10%, the value of $R_{0D}$ will decrease by approximately $1.478 \times 10^{-7}\%$. Conversely, if the parameter $\psi$ is decreased by 10%, the value of $R_{0D}$ will increase by approximately $1.478 \times 10^{-7}\%$.

The sensitivity analysis results, as depicted in Figure 3, reveals the sensitivity index of the parameter $\phi$ to $R_{0C}$, which is calculated at -0.62177. This means that if the parameter $\phi$ increases by 10%, the value of $R_{0C}$ will decrease by 6.2177%. Conversely, a 10% decrease in $\phi$ will result in a 6.2177% increase in the value of $R_{0C}$. Furthermore, the sensitivity analysis findings indicate a sensitivity index of $1.486\times 10^{-7}$ for the parameter $\phi$ to $R_{0D}$. In simpler terms, a 10% increase in $\phi$ will lead to $14.86\times 10^{-7} \%$ increase in the value of $R_{0D}$, while a 10% decrease in $\phi$ will correspondingly reduce the value of $R_{0D}$.

From the general sensitivity analysis, it becomes evident that the chemotherapy infusion rate ($\phi$) and the angiogenic dormancy rate ($\rho$) exert the most significant influence on the fundamental reproduction number $R_{0C}$. This implies that higher values of $\phi$ and $\rho$ will lead to a decrease in $R_{0C}$, thereby inhibiting the spread of tumor cells in gliomas. Similarly, among the nine parameters affecting $R_{0D}$, the chemotherapy infusion rate ($\phi$) and the angiogenic dormancy rate ($\rho$) are the most impactful. Notably, these impacts are in opposite directions. Increasing $\phi$ will result in the expansion of $R_{0D}$, leading to a wider tumor spread. Conversely, raising $\rho$ will cause a decrease in $R_{0D}$, thereby restricting the spread of glioma tumors.

4.   Results and discussion

In brief, this study unveiled a comprehensive model for glioma progression, emphasizing the efficacy of combination therapy and the role of antiangiogenic measures. Key parameters, particularly the chemotherapy infusion rate ($\phi$) and the angiogenic dormancy rate ($\rho$), significantly influence glioma proliferation. Sensitivity analysis identified $\phi$ and $\rho$ as crucial in decelerating glioma growth, shedding light on the interplay between drug-sensitive and drug-resistant cells. This insight was pivotal for refining treatment strategies and curbing disease progression. The significance of this study lies in optimizing therapeutic interventions through sensitivity analysis, providing valuable insights into glioma development and treatment effectiveness. Explaining these dynamics will contribute to advancing treatments for patients with glioma and propelling ongoing research in this field.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research is funded by the Directorate-General of Higher Education, Ministry of Education, Culture, Research and Technology, Republic of Indonesia under Master's Thesis Research Grant Number 3226/UN1 / DITLIT / Dig-Lit / PT.01.03/2023.

Conflict of interest

The authors declare there is no conflicts of interest.