A second dynamic population invasion study from reactive telegraph equation and boundary element formulation – a numerical assay about tumour growth in vitro

This paper is a continuation of a study already carried out on the use of the reactive-telegraph equation to analyse problems of population dynamics based on a formulation of the boundary element method (BEM). In this paper, the numerical model simulates the evolution of a tumour as a problem of population density of cancer cells from different reactive terms coupled to the reactive-telegraph equation to describe the growth and distribution of the population, similar to the two-dimensional in vitro tumour growth experiment. The mathematical model developed is called D-BEM, uses a time independent fundamental solution and the finite difference method is combined with BEM to approximate the time derivative terms and the Gaussian quadrature is used to calculate the domain integrals. The solution of the system nonlinear of equations is based on the Gaussian elimination method and the stability of the proposed formulation was verified. As the telegraph equation has a wave behaviour, a phase change phenomenon that can lead to the appearance of negative population density may occur, an algorithm was developed to guarantee the solution's positivity. Important results were obtained and demonstrate the effect of the delay parameter on tumour growth. In one of the tested cases, the results indicated an oscillatory behaviour in the tumour growth when the delay parameter assumed increasing values. The results of numerical simulations that sought to represent tumour growth, as well as the entire formulation of the boundary elements are presented below.

A second dynamic population invasion study from reactive telegraph equation and boundary element formulationa numerical assay about tumour growth in vitro Roberto Pettres a* , Andréia Assmann Pettres a , Eliandro Rodrigues Cirilo b a Federal University of Paraná, Brasil b State University of Londrina, Brasil * Autor correspondente (pettres@ufpr.br)

I N F O A B S T R A C T
Keyworks boundary element method reactive-telegraph equation with generation beginning and reproduction of cancer cells sufficient condition of the solution positivity tumour growth This paper is a continuation of a study already carried out on the use of the reactive-telegraph equation to analyse problems of population dynamics based on a formulation of the boundary element method (BEM). In this paper, the numerical model simulates the evolution of a tumour as a problem of population density of cancer cells from different reactive terms coupled to the reactive-telegraph equation to describe the growth and distribution of the population, similar to the two-dimensional in vitro tumour growth experiment. The mathematical model developed is called D-BEM, uses a time independent fundamental solution and the finite difference method is combined with BEM to approximate the time derivative terms and the Gaussian quadrature is used to calculate the domain integrals. The solution of the system nonlinear of equations is based on the Gaussian elimination method and the stability of the proposed formulation was verified. As the telegraph equation has a wave behaviour, a phase change phenomenon that can lead to the appearance of negative population density may occur, an algorithm was developed to guarantee the solution's positivity. Important results were obtained and demonstrate the effect of the delay parameter on tumour growth. In one of the tested cases, the results indicated an oscillatory behaviour in the tumour growth when the delay parameter assumed increasing values. The results of numerical simulations that sought to represent tumour growth, as well as the entire formulation of the boundary elements are presented below.

INTRODUCTION
In the last four decades, several therapeutic approaches to cancer treatment have been studied and the corresponding mathematical models of tumour growth have been developed (Byrne and Chaplain, 1995;Tao and Guo, 2006). Traditional tumour growth models are based on the diffusion equation (Preziosi, 2003;Swanson et al., 2003;Konukoglu et al., 2010;Painter and Hillen, 2013), as in heat diffusion problems, thermal conduction and dispersion of pollutants represented by the Fourier parabolic linear equation, which basically propagates information through space instantly with infinite velocity (Lurie and Belov, 2020). However, this equation has some shortcomings as it does not describe inertial effects and there are experimental evidences that the diffusive process takes place with finite velocity (Mittal and Dahiya, 2015). This behaviour is known as the 'Heat Conduction Paradox' and contradicts the so-called principle of causality, which states that information cannot travel faster than a finite velocity (Schwarzwälder, 2015).
In the case of a cancer evolution, using the diffusion equation to model the growth of a tumour is the equivalent of saying that every tissue or organ has cancer cells since the beginning of the oncological process, when in fact, the tumour can develop in nodules in a localized manner with limited morphology or in a progressive process. The tumour microenvironment is very complex in nature due to several simultaneous molecular mechanisms and simple diffusion can be a very simplified technique and may not represent the tumour microenvironment in detail (Sadhukhan and Basu, 2020).
L. Onsager in 1931, described that the Fourier's model was in contradiction with the principle of microscopic reversibility (Onsager, 1931), writing that this contradiction "is removed when we recognize that is only an approximate description of the process of conduction, neglecting the time needed for acceleration of the heat flow".
According to (Schwarzwälder, 2015), the Italian mathematician Carlo Cattaneo tried to overcome the problem of the infinite velocity of heat propagation by deriving a new equation to relate the heat flow Q and the temperature U, and therefore replace the Fourier law. Cattaneo introduced the characteristic thermal relaxation time , which is interpreted by Chandrasekharaiah as "the time interval required to establish constant heat conduction in a volume element, since a temperature gradient has been imposed on it" (Chandrasekharaiah, 1998)., that is, the time required to achieve thermodynamic stability. Therefore, the relaxation time introduces the idea of finite velocity of heat propagation, and it was first noticed by Maxwell (Maxwell, 1867).
Several attempts have been made to develop a precise mathematical model for heat, and perhaps two of the best known are the Maxwell-Cattaneo and Guyer-Krumhansl equations. The first introduces a relaxation time in the heat expression that has the effect of changing the governing equation to a form of wave equation, more preciously, the telegraph equation, which then exhibits a behaviour significantly different from the standard heat equation. Here, the most important related phenomenon is the so-called second sound, when temperature disturbance propagates like damped waves (Zhukovsky, 2016).
The second introduces non-local effects that incorporate interesting new phenomena, such as Cattaneo's Law of thermal viscosity and the hyperbolic heat equation. It is well adapted to the description of phonon gases where heat transport is not only governed by diffusion (like in Fourier's description) and second sound (like in Cattaneo's model in 1958) but in addition by ballistic transport, which is not the subject of this paper. For more details about Guyer-Krumhansl equation see (Zhukovsky, 2016).
In many references, these equations are also known as Maxwell-Cattaneo-Vernotte equations in honour of the French mathematician Pierre Vernotte, who derived the same equations almost at the same time (Vernotte, 1958).
In biological sciences, the telegraph equation can be used in studies of dispersal of species and/or populations in homogeneous environments in a probabilistic way, as in the work presented by Godoi (2021), by Malinzi (2021) about a model for oncolytic virus spread using the telegraph equation and also by Pettres (2021), in that the author present a variation of the telegraph equation, called the reaction-telegraph equation to simulate a problem of population invasion in a two-dimensional domain, making an analogy with the process of tumour growth in healthy tissue.
About tumour growth, the invasion model used to model the evolution of cancer indicates that once the tumour has grown to a certain size, it begins to actively invade healthy tissue (Hillen et al., 2015). An area of tissue can be considered as a 'cellular community' comprises 'species' of mesenchymal and epithelial cells in dynamic equilibrium with the environment and with each other (Gatenby and Vincent, 2003). A few tumour cells produced by successive non-lethal mutations begin to interact with the community of normal cells and are not recognized by this community, thus triggering the acquisition of space and the vital resources of the existing community, causing the growth of the tumour that may occur in two forms; active movement to nearby sites (invasion) or passive transport through the blood or lymphatic system (metastasis).
The tumour growth of few types of tumours are modelled by the exponential model (Johnson et al., 2019) given by equation (1), where u is the density of the cancer cells population and the parameter K1 defines the per capita growth rate or proliferation rate. Although exponential growth is a common initial assumption used to develop more complex models of tumour progression, few models strictly interested in characterizing tumour growth prescribe a fixed birth and death rate over time and population size (Johnson et al., 2019).
In the case of tumour dynamic saturation in the growth of various types of tumours is not well modelled by the exponential model. For this reason, this model applies only to avascular tumours, i.e., when angiogenesis has not occurred (Kerbel, 2000).
Indeed, tumour cells compete for oxygen and vital resources that is the reason why the logistic model presented by equation (2), tradtfits well in several cases ( 2) where the parameter K2 is the carrying capacity.
In the literature, one realizes that tumour growth does not follow a universal law (Varalta, 2014), however, according to (Gatenby and Vincent, 2003), two of the most used models are the generalized logistic model is given by, that describe well the growth of breast cancer (Spratt et al., 1996), and which the parameter K3 is the exponential decrease rate that represents an accelerated population degrowth when positive increasing values are used, and the Gompertz model given by, where K1, K2 and K3 are constants (variables u, K1, K2 and K3 are all positive).
The Gompertz model, while it does model the behaviour as a tumour increases in size, it is not an empirical model. The model has too many variables to consider, such as types of cancers as well as environmental conditions. These may even vary considerably for patients with the same types of cancers (Yorke et al., 1993). This model is the best to describe the volumetric growth in vivo (Michelson et al., 1987). For more details about mathematical models to study the various phases of cancer progression, see (Lowengrub et al., 2010).
One disadvantage of all the usual models is that they go to the carrying support faster than what is expected (Varalta, 2014). This is one of the advantages of the telegraphic model proposed by (Pettres, 2021), and it is here that this paper seeks to contribute with the theme, continuing the research based on a probabilistic model of an avascular-tumour growth in vitro simulating a tissue healthy (Ward and King, 1997), as well as a phenomenon of population invasion, under a more complete approach in relation to the model proposed by (Forys and Marciniak-Czochra, 2003), in which the authors used a simple logistic equation with time delay and simulated a diffusion problem.
For this purpose, an expanded version of the two-dimensional formulation of the boundary element called D-BEM presented by this author in (Pettres, 2021) is used, which is now presented as a continuation of the study already carried out on the telegraphic-diffusive phenomenon, here applied to tumour growth in vitro from different reactive terms defined by four growth models in the formulation.
D-BEM formulation used in this paper, widely deduced and explored in (Pettres, 2021), uses a time-independent fundamental solution to the potential problem in 2D combined with a timeadvance scheme based on the finite difference for the diffusion term (first derivative) and the wave term (second derivative) to solve the telegraph equation. This formulation reached a correlation level higher than 0.9 measured at the R 2 coefficient and with an average percentage error of approximately 8 per cent as described in (Pettres, 2021) from two problems with known analytical solutions. D-BEM formulation appears as a good option for the development of a hybrid numerical method, with the possibility of inserting different reactive terms in a simple way.
It is desired with this study, whose numerical results are reliable and can, indirectly, be used to evaluate the morphology of the tumour and the growth dynamics, explore an explicit threshold of the intensity of the immune response to control the tumour and to investigate the effects of the delay of reaction-telegraph equation in this new tumour growth model, with aims to contribute to the treatment of this complex disease, simulating from the first carcinogenic stimulus to the onset of neoplasia.

In vitro experiment and geometric modeltissue
The in vitro experiment is an important tool in cancer research, enabling the identification of carcinogens, the development of cancer therapies, drug screening, and providing insight into the molecular mechanisms of tumour growth and metastasis. Controlled simulations are described in (Ludwig, 1983), in which the authors use interferons to stimulate the tumour growth process from an in vitro assay.
The complexity of the model is largely dependent on the objectives. For example, preliminary screening of anticancer drugs can be performed in cell culture. Studies of invasion and motility of tumour cells can be performed with cells embedded in an ECM (Katt et al., 2016). For mode details about in vitro model see (Katt et al., 2016).
The growth takes place in vitro in order to eliminate the complexities due to vascularization and immunological responses, as noted in (Sutherland and Durand, 1973).
It is considered that there are nutrients such as sufficient oxygen and glucose available to all cells in the colony and possible declines in tumour growth will be considered as the result of chemicals produced by the eventual lysis of the cells due to the insertion of a specific treatment or drug from the parameter delay. Inhibition effects such as those that occur soon after the establishment of necrosis are neglected, as very little is known about this phenomenon, but Bertalanffy, cited by (Burkowski, 1977) estimates that, if such an effect is present, it is very likely insignificant compared to the inhibitory effect of chemicals from necrotic debris and therefore, it is chosen to ignore any inhibitory effects arising from viable cell metabolism. In the Invasion assays in vitro, a thin layer in thickness up to 1 mm of ECM is deposited on the porous membrane to model the basement membrane of the vasculature, typically in Matrigel although collagen and laminin are also used (Katt et al., 2016). For these reason, as an approximation, the geometric model used is a square tissue with two-dimensional domain = , ( , ),0 ≤ ≤ 1 , 0 ≤ ≤ 1 , illustrated in Figure 1, where the domain boundary is defined as  and it is subdivided in boundary elements n , simulating a two-dimensional in vitro experiment environment. Assuming that each cancer cell occupies an area of 1×10 -6 cm 2 , then the geometric model will contain a maximum of 1×10 6 cells.
The portion of cancer cells inserted in the experiment may vary according to the intended analysis. In the numerical tests that will be performed, cancer cells will be inserted into healthy tissue from a source geometrically represented by a Dirac delta function. If there is no growthinhibiting agent and there are also unlimited supplies, the transformation of the cancer will keep the number of cancer cells continuing to increase without limit. It is in this context that the telegraph equation, whose wavefront behaviour describes the population scattering or population invasion phenomenon, associated with a reactive term that contains in itself variables such as the growth rate or proliferation rate, carrying capacity and decrease rate can be useful to more realistically represent the growth dynamics of a tumour.

Reactive-telegraph equation
According to (Cirilo, 2019), equation (5) can be used to represent population evolution in a restricted domain, thus, the variable u represents population density and F(u) is a reactive term that describes the population growth. In this study F(u) is given by exponential model, logistic model, generalized logistic and Gompetz model presents in the equations (1), (2), (3) and (4), respectively.
The reactive-telegraph equation is deduced from the continuity equation and the constitutive equation of Cattaneo's (Méndez et al., 2010) in a probabilistic way (Godoi, 2021), taking the following form: In this paper, using a variation of the reactivetelegraph equation presented by (Pettres, 2021), in which the author represents the beginning and reproduction of cancer cells in healthy tissue including a punctual source term in equation (5), the study on the growth of the tumour from equation (6) and is called the equation of the reactive-telegraph equation with generation (RTEG), as follows is the diffusion coefficient that explains the tumour invasiveness or measure of how quickly the organism/cancer cells disperse, is the organism's finite velocity and is the organism's rate of changing direction (Holmes, 1993), ( , ) is the velocity cellular proliferation or time variation of population, is the relaxation time or delay time that produces a retard at solution, 2 ( , ) 2 is the term that allows to describe the proliferation effect in form of a travelling wave and Fs given by, where ( __ ) is a Dirac's delta, __ and are the localization and magnitude of source, respectively, represents the beginning and reproduction of the cancer cells in an empty domain, or healthy tissue, for example. Particularly, if = 0 it is found the diffusion equation with dissipative term and heat generation that was studied by (Pettres et al., 2015).
Assuming the essential and natural boundary conditions prescribed for equations (8) and (9), respectively: where boundary functions __ ( , ) and __ ( , ) in physical interpretation represent the population density and the invasion in the normal direction on the boundary , respectively. represents the part of the boundary of where the boundary condition is imposed, it is the part of the boundary of which is imposed the boundary condition q, where = . The initial condition at t = t0 given by (10).

D-BEM formulation
Considering the boundary value problem described by the equation (6) in 2D defined on domain Ω bounded by boundary . The integral equation of the BEM formulation for the RTEG equation can be written as follows, where ( ) is a geometric coefficient (12) at the collocation point , given by In this formulation of boundary elements, the fundamental solution of the Laplacian operator was chosen, and, being a solution of a static problem, the terms that depend on time in the integral equation for the RTEG equation are approximated using a time-advance scheme based on finite differences. This choice allows to insert different reactive terms into the formulation in a simple way, making D-BEM a hybrid numerical method. Thus, the time independent fundamental solution used in this D-BEM formulation is given by (Pettres, 2020), where r=X - is the distance between field and collocation points. The derivative of the fundamental solution with respect to the normal direction to the boundary is given by, where n is the outward direction normal to the boundary.
For simplicity, the time derivative presented in equation (11) is approximated by the backward finite difference formula (15) and (16) (Morton and Mayers, 1994), Assuming that the external source acts on the entire Ω domain and applying equation (17) to all boundary nodes and internal points, the following system of equations is obtained: where H bb and G bb are submatrices which result from the boundary integrals; the submatrices H db , G db , M db and M dd are matrices resulting from the domain integration, and are matrices resulting from the domain integration with derivatie of reactive term , F bd and F dd are matrices resulting from reactive term in the domain integration, and are matrices resulting from source term in the domain integration and I is an identity submatrix. The first superindex indicates the position of the collocation point and the second one the position of the field point, with b meaning boundary and d meaning domain. The subindex T+1 indicates the time ( +1) , the subindex T the time and the subindex indicates the time ( 1) . In this work a constant value is assumed for the time step , estimated according to (Wrobel, 1981) given by equation Erro! Fonte de referência não encontrada., where is the boundary element size.
By grouping similar terms of the system of equation Erro! Fonte de referência não encontrada., one obtains: Transferring the column coefficients of the matrices on the right-hand corresponding to the unknowns, to the left-hand side of the equation, one has: Equation (21) can be rearranged as equation (20) and solved with the same time marching scheme used previously.
After imposing the initial and boundary conditions, the system of equations (21) is solved and the unknown u and q values at the boundary nodes and u values at the internal points are determined at the time tT+1 . The u values are updated and the problem solution continues, recursively.
In general, after imposing the boundary conditions, one has.
where T  0 and: • +1 is the vector of unknown nodal values at time tT+1; • A is the coefficient matrix that contains terms relating to H, G and M; • +1 is a vector that represents contribution from the known values at the time tT+1, • the contribution from the previous time tT and • −1 represents the contribution from the past time tT-1.
In the following cases, the solution of (21) was found with the Gaussian elimination, whose solutions produced residues in the order of the machine's precision, in addition to being a wellposedness mathematical formulation. As for the convergence of the iterative system of equations, it is a sufficient condition of convergence that the matrix is strictly diagonal dominant, ensuring the convergence of the succession of the generated values to the exact solution of the system. For more details, see (Mahooti, 2020).
Comment on the singular integrationssame as in the classical 2D BEM formulation. For more details, see (Pettres et al., 2015). The domain integrals were calculated using the twodimensional Gauss quadrature method as described in (Pettres, 2021) and illustrate by Figure 2 that represent the location of a triangular cell and some Gaussian points and weights.    (Pettres, 2021).

Control of the phase change phenomenon -Sufficient condition of the solution positivityequation
The telegraph equation has a wave behaviour, there is a phenomenon of phase change that leads to the appearance of negative population density, which in this type of application has no physical significance. In (Cirilo et al., 2019) it was presented that the existence of negative solutions also depends on the parameters adopted, especially when using a constant value for the time delay parameter, with a high possibility of negative population density.
An alternative presented in (Alharbi and Petrovskii, 2018) and (Cirilo et al., 2019) to avoid this problem is the use of the so-called cut-off condition, establishing a null value for population density, whenever negative values for population density arise. Another alternative presented by (Oliveira, 2020) is to update the delay parameter whenever negative values for population density arise. This update takes into account the probability p=1-h, for h representing a one-dimensional segment, of an animal moving in a certain direction.
As already defined,  is the organism's rate of changing direction or rate of reversal of direction. This probability of reversing direction in a given period of time is known as a Poisson process with intensity  (Chiu et al., 2013).
A Poisson process meets the following criteria (in reality many phenomena modelled as Poisson processes don't meet these exactly): i) Events are independent of each other; ii) the occurrence of one event does not affect the probability another event will occur; iii) the average rate (events per time period) is constant and iv) two events cannot occur at the same time.
When memory effects are taken into account, successive movements of the dispersive population are not mutually independent, so that there is a correlation between successive steps (Fort and Méndez, 2002). According to (Tilles and Petrovskii, 2019), for that the results obtained with the reaction telegraph equation or the Cattaneoreaction system to be positive and for the boundary conditions to describe a physically sensitive situation and preserve positivity, the boundary conditions must be of the Neumann type or Robintype subject to trigonometric transformation, making the problem well posed. However, for an initial condition of the reaction-Cattaneo system given by a discontinuous function, the question of its positivity, strictly speaking, remains open and an analytical proof remains a challenge.
In these three alternatives, when using the socalled cut-off condition when there is a negative population density, in the first, when defining a null value for this population density, a fraction of the mass of the system and/or of the total population is disregarded. In the second alternative, by adopting a variable time delay parameter, the telegraphic phenomenon is considerably reduced and the diffusive phenomenon is accentuated and events are not mutually independent. In the third alternative, the guarantee of the positivity of the solution is made by changing the boundary conditions of the problem and also depends on the initial condition. Basically, the first alternative alters the mass of the system, the second alters the velocity of dispersion and the third alters the boundary conditions of the problem, violating the original problem in relation to the issue of biological mobility.
In this paper, it is proposed the use of the principle of mass conservation maintaining the velocity of dispersion without alters the boundary conditions of the problem, in which, occurring a negative population density at a certain instant in the interval of analysis, these values are taken in module and prorated among all other points of the domain with non-null population as a collective cell migration in tumours of highly differentiated tissues (e.g. lobular breast cancer and epithelial prostate cancer, where invasion by individual cancer cells is rarely detected) (Friedl and Wolf. 2003). Figure 5 illustrates how the total mass or population density of the system is conserved maintaining the dispersion velocity and boundary conditions. The graphs in the Figure 5 show some negative values on the blue line and the positivity adjustment on the red line, conserving the total number of the population even using a constant value for the relaxation time coefficient. During the iterative solution of the system of equations (21), at various times there is no need to use the positivity adjustment, as shown in Figure 6.

Solving D-BEM formulations for the twodimensional in vitro tumour model
All tumours follow a standard growth pattern, growing fastest in the beginning and eventually reaching a maximum size. A key component of any in vitro tumour model is a source of cancer cells, which in this paper are represented using boundary conditions q, as an invasion flow from adjacent structures and as a source term in the last example of each case. Four examples were studied for each of the reactive terms F1, F2, F3 and F4, representing different types of tumour growth, adopting = 0.0013 cm 2 /day, for diffusion coefficient that explains the tumour invasiveness or measure of how quickly the organism/cancer cells disperse, k1=0.012 cells/day for the tumour cell proliferation rate (Friedl and Wolf, 2003;Swanson, 1999) and relaxation time coefficient τ = 0, 5, 10, 20, 50 and 100 measured in days, representing some treatment or medication that slows down tumour growth, resembling transwell method (Katt et al., 2016), which are used to assay multiple parameters, such as the relative invasiveness of different cells and the effect of drugs or gene manipulation on motility.
The first case to be analysed is the following: Case 1 -Considering the exponential growth for the reactive (reproductive) process, has The solutions of the equations systems were based on the Gaussian elimination method, whose results for D-BEM showed stability. Results of the formulation D-BEM for reactive-telegraph equation are illustrated below. Example 1.1 -Beginning and reproduction of cancer cells in the centre of healthy tissue in which cells do not invade neighbouring structures for the exponetial growth.
In the first example, the punctual source term in the domain is considered constant along time mensure in , given by: The value adopted for this production function is related to the maximum capacity of occupation of cells in the area defined in the experiment, 10 6 per cm 2 , with = 0 for an analysis period of up to 100 days.
The use of this function is strongly related to the existence of a problem in the cellular apoptosis system, which is the programmed death of the cellular connection. Normally, a cell that has passed through its life cycle dies. In its place, a new subpopulation of the cell cycle develops over time. But with the transformation of cancer, this natural mechanism is disrupted, as a result of which that cell does not die, but continues to grow and function in the body.
It is this internal mechanism that is the basic basis of tumour formation, which has a tendency to uncontrolled and unlimited growth. That is, in fact, this type of cellular structure is a cell that is not capable of death and has unlimited growth, assuming that the cells don't change their proliferation kinetics when implanted in the group of cells of the same cell line.
To solve the RTEG model, proper initial and boundary conditions are required. The initial condition is given by equation (26) representing null population density and null population variation only at the beginning of the reaction in the domain and the micro-environment is full of nutrients and the necrotic cells are not found in the domain of interest initially. representing null mobility at the border, indicating that the population generated by remains in the domain and also does not allow populations who may come from the vicinity of the domain to invade the tissue. This case represents a benign tumour, in which this tumour kind is well-defined local masses, whose cells do not invade neighbouring structures or spread at a distance through the blood and/or lymphatic stream, so that they do not form metastases.
The D-BEM results for this example are illustrated below. The source term leads to a considerable increase in population density in a short period of time and has unlimited growth due to the capacity of replication in the proposed system and the absence of the delay factor (τ = 0) in the analysis, as illustrated in the first image of the sequence in Figure 7.
In cases with ≠ 0, the concentration of cells started to show decreasing values when using increasing values for the relaxation coefficient t. It is important to highlight the effect of the delay in the concentration of cells observed in the graph to the right of each test in this example.
With the use of increasing values for the τ parameter, it is possible to note that there is a delay in the population growth and an oscillatory behaviour, as illustrated in Figure 7 with τ equal to 5, 10, 20, 50 and 100. From the theory of delay differential equations, (Forys and Marciniak-Czochra, 2003) state that for one equation with delay, an oscillatory behaviour is possible even if for the same equation without delay there are no oscillations, as in the purely diffusive case. This is a very important result, as it illustrates two phenomena, one of which is the origin of cancer cells that are constantly produced in the tissue according to an exponential growth and the effect of the delay parameter in expansion of number of these tumour cells. Depending on the type of cancerous evolution, a certain delay factor can be dimensioned so that there is enough time for the tissue to regenerate or even, to keep the tumour at a stable size, through continuous medication. Typically, very aggressive tumours take an average of 60 days to double your size and 100 days for nonaggressive tumours.
Analysing the concentration of cancer cells per cm 2 for each τ at certain times in time up to the limit of 100 days, the following table is presented.  Table 1 shows the rapid growth in the value of the concentration of cancer cells per cm2 in the reaction-diffusive case (τ = 0), with a homogeneous distribution through out the tissue in the interval of up to 100 days. With the use of increasing values for the delay parameter, there is a decrease in the rate of growth of the concentration of cancer cells, as illustrated in Figure 8, when we compare the cases with τ equal to 0, 5, 10, 20, 50 and 100, respectively.  In Figure 9, when using τ = 0 all points u(x, y, t) in the domain, they are changed instantly (diffusion equation -infinite velocity), however, when using τ ≠ 0, in this case, τ = 5 and τ = 10, respectively, the points in the domain furthest from the source of cancer cells gradually perceive their effect as a moving with delay in space over time and the density of infected cells increases in the radial direction due to the symmetrical geometric model, as illustrated by Figure 10. The results obtained with the RTEG model are similar to the in vitro experimental results performed by (Jiang et al., 2014) illustrated by Figure 11, whose cell culture invasion assay provided a physiological approach to assess tumour invasion and offered a visual component that can be quantified through image analysis.  In Figure 12, a geometric analysis was performed, defining four circumferences with radius r1, r2, r3 and r4, repectively, which were used to compare the space occupied by the population of cancer cells in the numerical model and in the in vitro experiment. It is possible to observe that the numerical results provide a good approximation of population growth and spread, similar to the real experiment. Note that with the use of the telegraph equation, the tumour growth process can be significantly represented and that the delay coefficient can be scaled depending on the type of cancer. At this moment, it is possible to afirmate that, biologically, the delay coefficient has effects on the tumour growth and on the spatio-temporal distributions of tumour cells and telegraph model to tumour growth shows is promising.
Deepening the analysis of Figure 7, an important and unexpected behaviour in tumour growth is revealed after the fiftieth day, approximately, in which a great growth is observed in all cases with τ ≠ 0 , predicting that the tumour growth may significantly accelerate its process of cell invasion. Table 2 shows how the number of cancer cells grows from the previous period. In Table 2, despite the density of cancer cells growing up to 8 times between 5 and 10 days, in this period there is still a few invaded cells and any growth represents a rapid growth, however, This modelfor population growth was first proposed by Verhulst in 1838 (Bacaer, 2011). The dynamics are modelled by the first order ODE, wherethe parameter k1 = 0.012 cells/day defines the per capita growth rate or proliferation rate and parameter k2 = 10 6 cells/cm 2 is the carrying capacity defined for the in vitro model based on the diffusion model obtained for 100 days in the case 1. For u ≪ K2 the growth rate is K1, but as u increases a quadratic nonlinearity kicks in and the rate vanishes for u = K2 and is negative for u > K2. The nonlinearity models the effects of competition between the organisms for food, shelter, or other resources. There are two fixed points, one at u = 0, which is unstable f′(0) = K1 > 0.
The derivative of the reactive term in this example is: From equation (17) Note that in equation (29) there are terms of the form 2 ( , ) , which in this study was chosen to start them in the iterative process at t and − , as variables known in previous times. The matrix system of equations is similar to equation (21). In the next step, the study for the logistic model as a reactive term begins. Example 2.1 -Beginning and reproduction of cancer cells in the centre of healthy tissue in which cells do not invade neighbouring structures for the logistic growth In the second example, the punctual source term in the domain Ω is defined as equation (25), initial conditions (26), boundary conditions (27) and τ equal to 0, 5, 10, 20, 50,and 100. The D-BEM results for this example are illustrated below. In the first analysis with τ = 0, it is possible to verify that there is no wave behaviour, only diffusion with reation. In analyses performed with τ ≠ 0 illustrated by Figure 13, results were obtained similar to the results obtained with the exponential reactive term explored in Case 1 (Figure 7). Figure   14 compares the D-BEM results for Case 1 and Case 2 with τ = 50 and K2 =10 6 cells/cm 2 , where it is possible to verify the similarity between the two tested cases, for example. In a complementary analysis presented in Table 3, values 100 times lower than K2 and 100 times higher than K2 were used for the carrying capacity coefficient. Comparing the values in Table 1 (τ = 50 and K2 =10 6 cells/cm 2 ) with the values in Table 3 (K2 =10 6 cells/cm 2 ), it can be seen that the variation between the population density values after 10 days is approximately 1%, indicating that, although different coefficients that represent the carrying capacity present in the logistic model that is used as a reactive term, the results are similar to the case with an exponential reactive term with τ equal to 50. These results indicate that the resulting estimation of the carrying capacity K2 is biologically irrelevant, and it can be reaffirmed that, based on the simulated conditions, the logistic model resembles the exponential growth model.
Case 3 -Considering the generalized logistic model for the reactive (reproductive) process, has According to (Aviv-Sharon and Aharoni, 2020), a little more than half a century ago, an extension of the classic logistic model, allowing for more flexible curvature of the S shape where the growth curve is asymmetrical, was introduced, establishing the so-called Richards curve or generalized logistic model proposed by Richards in 1959(Richards, 1959. Inspired by population biology, this model assumes an initial phase of exponential growth that saturates as the number of cells reaches the capacity of space and available nutrients.
Again adopting K1 = 0.012 cells/day for the tumour cell proliferation rate, K2 = 10 6 cells/cm 2 for the carrying capacity, relaxation time coefficient τ = 0, 5, 10, 20, 50 and 100 measured in days and K3 = 2.5 for the exponential decrease rate. The derivative of the reactive term in this example is: From equation (17) Note that in equation (32) there are terms of the variable u in which, when the multiplication present in some domain integrals is performed, they result in terms of the shape u K 3 (X,t) and u K 3 +1 (X,t), which leads to this stage of the study considering a population growth based on fractional telegraphic reaction, which in this study was chosen to start them in the iterative process at and − , as variables known in previous times. Equation (32) can be reorganizedin a matrix system of equations similar to equation (21), the results of which are presented in example 3.1.
Example 3.1 -Beginning and reproduction of cancer cells in the center of healthy tissue in which cells do not invade neighboring structures for the generalized logistic growth In the third example, the punctual source term in the domain Ω is defined as equation (25), initial conditions (26), boundary conditions (27) and τ equal to 0, 5, 10, 20, 50 and 100. The D-BEM results for this example are illustrated below. In the first analysis with τ = 0, only the diffusive behaviour with generation in population growth is perceived. However, in the rest of the analyses performed with τ ≠ 0, like illustrated by Figure 15, results were obtained similar to the results obtained with the exponential reactive term explored in Case 1 (Figure 7). In a complementary analysis presented in Table 4, was used K3= 1, 1.5, 2, 2.5, 4, 5.5 and 7, for the exponential decrease rate. Comparing the values in Table 4 to different K3, it can be seen that the maximum variation between the population density values is approximately 0.03% (for K3 = 1 and 7 in 100 days), indicating that, although different coefficients are used to represent the exponential decrease rate, present in the generalized logistic model, the results are similar to the case with an exponential reactive term with τ equal to 50. These results indicate that the resulting estimation of the exponential decrease rate K3, when K3 is between 1 and 7, is practically irrelevant for a period of time up to 100 days.
The Gompertz model is characterized by an exponential decrease of the specific growth rate with the rate denoted here by K1, which is a constant related to the proliferative ability of the cells. K2 is the carrying capacity.
According to (Vaghi, 2019), the tumour growth is not entirely exponential, provided it is observed over a long enough time frame (100 to 1000 folds of increase) and the specific growth rate slows down and this deceleration can be particularly well captured by the Gompertz model. Despite its abillity, the initial rate of tumour growth is difficult to predetermine given the varying microcosms present with a patient, or varying environmental factors in the case of population biology. In cancer patients, factors such as age, diet, ethnicity, genetic pre-dispositions, metabolism, lifestyle and origin of metastasis play a role in determining the tumour growth rate. The carrying capacity is also expected to change based on these factors, and so describing such phenomena is difficult. However, in this paper K1 and K2 are defined as constant, that is, k1 = 0.012 cells/day and k2 = 10 6 cells/cm 2 , in order to compare the results with the values obtained in cases 1, 2 and 3, because during the preliminary tests it was found that the growth based on the Gompertz model is highly sensitive to variations in the carrying capacity coefficient.
The derivative of the reactive term in this example is: From equation (17)   In the fourth and last Case, the punctual source term in the domain Ω is defined as equation (25), initial conditions (26), boundary conditions (27) and τ equal to 0, 5, 10, 20, 50 and 100. The D-BEM results for this example are illustrated below. In the first analysis shown in Figure 16 with τ = 0, the population growth behavior resembles an exponential curve, reaching a slightly smaller amount of cancer cells, when compared to the cases with the reactive terms F1, F2 and F3 with null τ.
In tests with τ different from zero, a logarithmic behaviour is perceived, differently from the behaviour presented for the other three cases. In this analysis, an increase in the number of cancer cells can be seen similar to an S-shaped curved growth (sigmoid curve), common in logistic models that take into account the carrying capacity K2.
Furthermore, in the rest of the analyses performed with, as shown in Figure 14, results were obtained that are not similar to the results obtained with the reactive terms F1, F2 and F3, reaching lower values for the same time interval.
In two complementary analyses presented in Table 5 and 6, values 1000 times lower than K2 and 1000 times higher than K2 were used for the carrying capacity coefficient assuming τ = 50 and in the second analysis, the results shown in Figure  16 are presented in detail. Using the reactive term based on the Gompertz growth model, the effect of the carrying capacity is perceived as shown in Table 5, in which, the growth is highly influenced by the carrying capacity, growing at a higher rate and reaching very high values after 50 days approximately. Table 6 presents the cancer cell counts from the tests performed in this section. Results with the parameter τ different from zero show a slow growth until the 25th day, however, after the 50th day, there is a rapid increase in cell count, but on the 100th day, they present values lower than Cases 1, 2 and 3. This result is explained by the effective effect of the carrying capacity parameter. Important information contained in the graphs in Figure 16 is that, as the growth curve occurs in a logarithmic way, the population density tends to remain constant after a certain period of time, reaching a plateau. This affirmation was proven by simulating the same problem over a longer period of time, 365 days, as shown in the Figure 17  In Figure 17, the effect of the delay parameter τ is evident, causing the cancer cell population count to oscillate over time, converging to a limited amount, unlike the case with τ = 0 in which growth occurs almost exponentially reaching values of the order of 10 7 . In this analysis, convergence to a limited value is observed with a time interval less than 150 days for cases with τ = 5 and τ = 10. For cases with τ = 20 and τ = 50, convergence was observed after approximately 200 days and for the last case, with τ = 100, the amount of cancer cells stabilized after 300 days. The red dashed line is used to highlight the value to which the population converges, reaching lower population values when increasing values are used for the parameter τ.
Using the reactive term based on the Gompertz growth model, it is observed that the growth is smooth at the beginning of the reaction, reaches a maximum value, starts to decrease due to the carrying capacity, returns with growth and then stabilizes forming a plateau. This is an important result and differs from the three cases tested previously, because in cases 1, 2 and 3, what was noticed was that the telegraph equation delayed population growth and here, with the Gompertz model, in addition to the delay, also population regrowth, resembling a real organism undergoing treatment, in which growth is inhibited by some physical, chemical or biological agent (radiotherapy, chemotherapy or hormonal therapy). This agent may be related to the delay parameter τ, in addition to being useful in tissue recovery analysis, which would give enough time for new cells to be generated and replace invaded tissues, for example, or lead to surgical eradication, without causing damage to surrounding tissues in the case of benign tumors and the use of radiotherapy, chemotherapy and hormone therapy in malignant tumors, which can be used to alleviate signs and symptoms and try to prevent cancer recurrence.

General comments
To reinforce the claims made in this paper about the use of RTEG to model the dynamics of a tumour in vitro, in (Jiang et al., 2014) the authors have identified different migratory patterns from the clinical images of the tumours, which were collected in Sun Yat-sen University Cancer Center. They have observed that the invasive liver tumour spread in faster rate (characteristic of superdiffusion or even ballistic diffusion) than the adrenal tumour (characteristic of sub-diffusion). From this phenomenon, they have concluded that the rate of tumour expansion is varied rapidly depending upon its niche. So, the simple diffusion is not suitable to model these types of phenomena as it converges quickly, and this way, is here that which this paper can contribute, when considering the use of the telegraph model in tumour growth.
If the reactive-telegraph equation with generation were used to model the evolution of cancer cell reproduction, for example, the use of the reactive terms F1, F2 and F3 show similar results and, in a way, would be it is convenient to use the exponential growth model (F1), as it is the easiest to implement computationally. However, the tests carried out with the Gompertz model presented possibilities not observed in the other growth models studied in this work, highlighting the possibility to simulate a fast growth, a decrease, and a stabilization of the population density of cancer cells. This decrease and stabilization of the amount of cancer cells can be modelled with the delay parameter τ, which would act as a growthinhibitory agent, taking into account some type of specific treatment or medication, which, in addition to being useful in the analysis of tissue recovery, would also provide enough time for new cells to be generated and to replace the invaded tissues, for example, avoiding the failure of some organ or tissue.
The results of these four cases reveal the importance of this numerical study, which, by means of mathematical modelling, can assist in the understanding of cancer staging, which is very important because it determines the most appropriate treatment and is an important indicator of prognosis and possible complications, in addition to assessing the cellular physiology underlying a tumour and predicting invasive behaviours in healthy tissues and in assessing the response to a treatment, just by changing parameters in the numerical model.
Regarding the condition of the solution positivity, it was not necessary to use this condition in the study with the reactive term based on the Gompertz model, revealing a curious fact regarding the wavefront behavior of the invasion process, indicating that there was no reflection of the waves at the domain boundary, evidenced by the population accumulation in this region provided by the Neumann-type boundary condition, which in the context of this application, allows tumor growth at a variable and decreasing rate, limited by the coefficient of carrying capacity.

CONCLUSIONS
In this article, the choice of the fundamental solution of the Laplacian operator (eq. 13) and the use of finite differences (eq. 15 and 16) proved to be a good option for the development of a hybrid method to solve the non-homogeneous telegraph equation, whose insertion of different reactive terms was made in a simple way in this mathematical formulation.
Using the reactive-telegraph equation and of the D-BEM formulation, it was sought to numerically simulate the population growth dynamics of a group of cancer cells from a two-dimensional model that resembles a controlled assay of tumour growth in vitro.
Four population growth models were used for the term reactive, in which, the exponential model, logistic and generalized logistic models, presented similar results. In these three cases, it was noted that as increasing values are used for the relaxation time coefficient, the effect of the second sound is perceived and the replication of the population's cells becomes slower. This damping feature can be dimensioned to consistently represent the growth of some types of cancer and, mainly, to relate it to a specific type of treatment or medication, which causes a delay in tumour growth.
The results obtained with the fourth model tested, the Gompertz model, indicated that in addition to a delay in population growth, a decrease and stabilization of the population was also observed, strongly influenced by the use of the τ parameter and the coefficient of carrying capacity.
The models developed allowed simulating different tumour growth scenarios and different delay parameters were adopted. This model was investigated from a phenomenological point of view, measuring some parameters to characterize the growth of the avascular tumour over time. The tumour micro-environment is very complex in nature due to several competing molecular mechanisms. Diffusion can be an overly simplified technique and does not reflect the detailed view of the tumour micro-environment when some type of treatment is inserted. In addition, as the simulation parameters can be modified due to different biochemical and/or biophysical processes, the model presents sensitively the variations in the results, guaranteeing the robustness of the model.
As desired with this study, based on the developed mathematical formulation and the numerical results obtained, it is possible to evaluate the tumour morphology and the growth dynamics, allowing to explore the effect of the variation of parameters related to tumour growth, as the parameter K1 that defines the per capita growth rate or proliferation rate and K2 that defines the carrying capacity, present in the reactive terms of the equation telegraph, which may be related to a particular type of treatment. In this study, the parameter K3, which defines the exponential decay rate, was not relevant in the population count.
This can help us understand how an avascular tumour interacts with its micro-environment and what kind of changes in its growth can be seen if a drug or clinical treatment is administered before angiogenesis.
It is recognized that, in order to understand the evolution of cancer, many still needs to be studied on first, but this paper will likely help to develop better mathematical models to simulate diagnoses, therapies and prognoses, as this formulation allows to simulate a real problem, taking into account the adjustment of the growth parameters involved, the control of cell proliferation and migration mechanisms using specific boundary conditions and by the frequency of growth stimulation of a certain type of tumour cell by adjusting the Fs source.
Finally, it is emphasized that the models analysed are not restricted to a specific system and can thus be useful in the study of wave fronts in biophysical and chemical problems where there is the presence of reaction-diffusion, as in nerve conduction, in the growth of colonies of bacterial, population dispersion and tumour growth models in healthy tissue, as widely discussed in this paper, especially considering the importance of latent time for most pathologies. For the future, based on the results of this second study, it is intended to direct the research towards comparative applications taking into account real data of tumour growth in vitro and the results of the numerical model developed to model the parameter τ as a growthinhibiting agent, taking into account some type of treatment or specific medication.