Implications of the “Subquantum Level” in Carcinogenesis and Tumor Progression via Scale Relativity Theory

1. Introduction

The last 25 years witnessed tremendous achievements in cancer diagnose and treatment. Technology currently permits small size tumors (like breast cancers) diagnosis and treatment, ductal cancer in situ currently including 25...30% of all freshly diagnosed breast cancers at the majority of medical centers [1]. Thus, early detection now allows the understanding of growth patterns. Surgeons are in the front line of technological and basic scientific medical advances. Current ideas, such as the physiological characteristics of shock, organ transplantation, antisepsis, wound healing, or sequence medical care, are cast by surgical investigators.

The field of mathematics suffered an identical evolution. Revolutionary mathematical branches such as topology, fractals, chaos theory, and development of nonlinear descriptive strategies have provided mathematicians new inventive tools to create growth models and to behaviors at the small environmental level [2, 3]. Growth, angiogenesis [4], cell-to-cell adhesion [5], hydrogen ion concentration regulation and drug delivery [6] can now systematically be described using specific formulas. From a clinical viewpoint several of these formulas could seem simple, however they put together a very important foundation for descriptive insight.

What is currently lacking is a connection between these two naturally and mutual analysis endeavors. For oncology surgeons, the ability to mathematically analyze and predict patterns of growth provides precise techniques that are beneficent for both current and future therapies. For mathematicians, defining the clinical factors essential for growth development and metastasis can provide realistic insight into these biological processes, successively allowing the event of correct, clinically relevant mathematical formulas. In almost every dedicated medical institutions, the teams that have comprehensive cancer are being led by surgeons. Mathematics that can be applied in the oncology field provide a chance to expand the leadership role of the surgeons and to raise awareness about the significance of understanding growth behavior and, also improve cancer treatments.

The present study aims at defining a new concept of carcinogenesis and tumor progression. Consequently, we use the natural ‘environment’ where malignant tumors grow, space(-time) with non-integer fractal dimension, questing for further applications of the newly discovered and intriguing phenomenon of tumor self-seeding by circulating cancer cells (CTC). More precisely, we assume that the metastatic tumor cells move (through the systemic circulation, yet not necessarily only there) as a coherent wave, or even more precisely, a chemically pumped travelling laser wave with oxygen. The extracellular matrix (ECM) and in particular, the tumor microenvironment (TME) are assumed as non-differential media endowed with holographic properties and may be good candidates for “recording” materials. As a result, the tumor self-seeding by CTC may be proved mathematically, the fact that the CTC returning to the initial tumor site and fueling the primary tumor growth or even grow a new tumor is a particular case of complete holography (i.e. a hologram which does not represent only the virtual object’s image, but it becomes the very object - which we believe, is a characteristic of the living organisms). We believe our findings may provide new opportunities to set up new targeted therapies that may slow down or even prevent tumor progression.

According to the Pribram-Bohm’s holographic theory (http://en.wikipedia.org/wiki/Holonomic_brain_theory) of the brain, the intercellular and intracellular communication implies the existence of a fractal medium equivalent to the vacuum between the elementary particles. Since vacuum dynamics is studied using quantum mechanics, it is only natural that the status of the fractal medium implies the use of a quantum type mathematical formalism (i.e. Scale Relativity theory either in its Schrödinger type representation or the fractal hydrodynamic one) applicable to different resolution scales (mesoscopic for intercellular communications or nanoscopic for intracellular ones). In Nottale’s interpretation, each resolution scale is characterized by a Planck type constant $$\overline{h}$$

$$\overline{h}=\hslash {\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}$$

where ħ is the standard Planck’s constant, dt is time’s resolution scale, τ the time’s reference scale and D_F is the fractal dimension of the motion curve.

2. The biology of cancer

2.1. Cancer, what should be noticed

Cancer or malignant neoplasm is a class of diseases that rises from the anomalous behavior of normal tissue. Cancer cells are aberrant cells which have acquired malignant traits such as uncontrolled growth (cells continuously proliferate), tissue invasion (they intrude into normal tissue and destroy it) and metastasis (they spread outside the location of the body where they were originally generated). Additionally, the term tumor or neoplasm is used to indicate an abnormal swelling of tissue caused by an excessive cell proliferation.

A tumor can be of benign or malignant nature, while benign tumors are self-limiting, do not express patterns of invasion, and they do not metastasize, malignant tumors do possess all these characteristics. The term malignant tumor is also used as synonym for cancer, although some cancers, such as leukemia, do not form tumors.

Cancer cells develop these malignant features because of genetic mutations, accumulated during the organism lifetime. Cancer is in fact a multi-step chance process that transforms a normal cell into a tumor cell, after having collected a set of 5...8 crucial genetic alterations [7-9] as schematically shown in Fig. 1.

A newborn malignant cell, expressing aberrant traits, can lead to the formation of cancer and, in most of the cases, of a tumor. Without treatment, the destructive behavior of such colony of cells is usually lethal for the patient. The probabilistic nature of this disease and the increase in life expectancy had made cancer the second cause of death in the industrialized countries (see any cancer statistics). Nevertheless, cancer it is not a modern disease and it was known since the antiquity: Egyptians of the New Kingdom [10], Greeks [11] and Romans [12] accurately described medical treatments for tumor removal. It is only within the last two centuries however, that due to the higher standards of living, cancer has become one of the main life-threatening diseases.

Figure 1.

Acquisition of the tumorigenic phenotype by a population of normal cells through multiple genetic mutations.

2.2. Distinguishing traits of cancer

The tumorigenic properties, generically discussed in the previous section, have shown to be common to almost all cancers. They have been studied since the dawn of cancer research and they can be enumerated and defined with a relatively high accuracy. These hallmarks are a set of characteristic traits typical of cancer cells that are essential for the formation of a macroscopic malignant neoplasm [13]:

Self-Sufficiency in Growth Signals. All cells communicate through signals. A biological signal is, in most of the cases, a protein able to deliver a particular piece of information by binding uniquely to specific receptors on the cell surface. Normal cells need mitogenic growth signals to proliferate (signals that allow and stimulate cell proliferation). Those signals are regulated by the homeostasis of the tissue and they guarantee a correct balance between cell proliferation and death, according to the needs of the organism. In order to lead to cancer, tumor cells may develop the ability of self-generating such signals in one way or another. One possible way is a genetic aberration in one of the fundamental genes responsible for the building of the signaling pathway, for instance the RAS oncogene [9,14]. As consequence, the associated component of the signaling system would become constitutively active and hence, independent by the signal molecule. A second option is the self-production of growth factors that would stimulate growth by paracrine signaling, where a cells stimulates the neighbors and vice-versa or even autocrine signaling, when the cell stimulates its own receptors as shown in Fig. 2.

Figure 2.

Example of self-signaling (autocrine): the cell produces its own growth factors which stimulate the growth receptors on the surface (Dr. W.H. Moolenaar, Netherlands Cancer Institute).

Insensitivity to Antigrowth Signals. As counterparts of growth factors, homeostasis employs growth inhibiting signals as well. These signals act similarly to their antagonists but they promote cell cycle arrest or cell quiescence, rather than proliferation. An example of a crucial gene involved in anti-growth pathways is the retinoblastoma protein (pRb). The retinoblastoma protein is capable of altering the function of the E2F transcription factors and control the expression of the bank of genes essential for the transition from GAP-1 phase to DNA Synthesis phase of the cell cycle [15]. The disruption of such pathway results in the insensitivity of the cell to anti-growth signals.

Evading Apoptosis. Apoptosis is a mechanism of controlled cell death. Through special signals, a cell has the capacity of terminating itself in a highly regulated way. A normal cell dying by apoptosis undergoes a sequence of events such as condensation, fragmentation and phagocytosis. This avoids the cell to free potentially dangerous enzymes and proteins stored inside its cytoplasm and its nucleus. During apoptosis the cell membrane is kept intact while, in 30...120 minutes the cell is fragmented in small parts or apoptotic bodies, still protected by pieces of membrane. Those cell leftovers are successively phagocytated by macrophages within the next 24 hours [16]. Apoptosis is a common mechanism of cell death and takes part in the homeostasis of a healthy organism as well as in its embryogenesis and in its morphogenesis. When any cell violates such homeostasis, an apoptotic signal is delivered to it. Therefore, in order for cancer cells to develop into a malignant lesion, it is necessary to deactivate apoptotic signal pathways. A mutation in the p53 tumor suppressor gene (TSG) is one of the most common ways to acquire resistance to apoptosis because p53 regulates the whole signaling process of programmed cell death. Indeed, more than 50% of human cancers carry a mutation in the p53 tumor suppressor gene [17].

Limitless Replicative Potential. Even with all the anti-growth and anti-apoptosis pathways triggered off, a cell could not generate a vast population able to form a tumor. That is because of the intrinsic proliferation limit of all mammalian cells. All chromosomes have an ending cap called telomere, a T-loop non-coding DNA sequence (2...50 Kb) that prevents the end of the chromosomes from attaching to other genetic material. At every mitosis the cell loses a small part of its telomeres because of the impossibility for DNA duplication enzymes, for instance DNA-polymerase, to continue working until the very end of the genome (Fig. 3). This limitation is due to the fact that enzymes like DNA-polymerase always move in the 5’...3’ direction of the DNA sequence, so when the side of the replication is opposite, a small part of the genome is lost. The shortening of the telomeres induces cell senescence, a state of cellular elderly where division no longer occurs. This avoids genetically unstable cells to replicate. Senescence starts after the so called Hayflick limit [18] of about 50 cell divisions. In cancer cells instead, the disabling of the pRb and the p53 pathways allows unlimited replication, until the point when the telomeres are completely absent. Once having entirely consumed the telomeres, the cell population is believed to undergo a phase of massive genomic instability, causing extended cell death. The high selective pressure induced by this crisis may permit specific resistant clones to emerge (Fig. 4). Those survivor cells would be immortalized (unlimited proliferative potential) by finding ways to maintain their telomeres long enough. A possible way is the over expression of the telomerase gene [19] which appears to take place in 85...90 % of cancers. Telomerase is a telomere-rebuilding enzyme normally expressed in germ line cells and stem cells, in which immortalization is an essential feature. Once immortalized, malignant cells have made a further step towards the formation of cancer.

Figure 3.

Illustration of the end replication problem: at both sides of the copying, the leading DNA strand has lost part of the telomeric sequence, which stops at the 5’ end of the parental strand, whereas the lagging strand results completed until the very end (Dr. R. Beijersbergen, Netherlands Cancer Institute).

Sustained Angiogenesis. In order for a cell to survive normally, it must rely within 100 μm from a capillary blood vessel [13]. For this reason, the initial exponential growth of a newborn malignant neoplasm causes a shortage of nutrients among cancer cells. Local pre-existent vascularisation is never enough to sustain growth for more than 10⁸ cells. The colony must therefore develop angiogenesis-triggering capabilities [20, 21]. Angiogenesis is the process of formation of new blood vessels in response to a stimulus secreted by poor vascularised tissues. Angiogenesis is important for the organism morphogenesis and even later maintains the correct supply of nutrients for all tissues. Fast growing cells, such as cancer cells, start soon to starve, and have the need of additional blood supply in order to keep expanding. A possible solution is the production by cancer cells of vascular endothelial growth factors (VEGF) and fibroblast growth factors (FGF1/2) which bind to the transmembrane receptors of endothelial cells (cells covering the interior surface of blood vessels) stimulating their growth towards the signal concentration gradient [22]. Angiogenesis is the principal mechanism that transforms a microscopic malignancy into a macroscopic tumor and, also in the later stages, it is necessary for a lesion to grow and sustain itself. This implies that angiogenesis is an important target for anti-cancer drugs like thrombospondin-1 [23] and bevacizumab [24], also known as avastin.

Figure 4.

The progressive shortening of the telomeres leads to a massive cell death due to the induced genomic instability (death by genomic catastrophe while duplicating). From such process of intense genetic mutation and selection an immortalized clone could emerge (Dr. R. Beijersbergen, Netherlands Cancer Institute).

Tissue Invasion and Metastasis. The most dangerous and destructive features of cancer are tissue invasion and the consequent metastasis. Its ability of forming distant colonies or metastases all over the body represents the cause of 90% of all cancer related deaths [25]. Normal cells are usually unable to travel outside their own tissue due to their necessity to be anchored and reside among similar cells. An eventual detachment from the extracellular matrix or ECM (a complex structure of proteins and specific cells forming the tissue scaffold and microenvironment – see Sec. 6.1) would occur in a form of apoptosis called anoikis [26]. Contrary to their normal counterparts, cancer cells are able to survive the loss of anchorage, to travel through the vascular system and form distant tumors elsewhere (Fig. 5). The traits expressed by invasive and metastatic cancer cells are principally loss of cell-to-cell adhesion, anchorage-independence, chemotaxis (migration towards a diffusible substance gradient), haptotaxis (migration towards a non-diffusible substance gradient) and production of matrix degrading enzymes (e.g. Matrix metalloproteinase) which cleave the extracellular matrix [27-29] making space for invasion and freeing growth and angiogenic factors trapped inside.

Figure 5.

Tissue invasion is a multi-step process that requires the cancer cell to have developed many malignant traits, necessary for the formation of new distant colonies called metastases [27].

3. Mathematics of cancer

In comparison to biology, cell biology, and drug delivery analysis, mathematics has, to date, made comparatively very few contributions to this field of research. A statistical analysis of the PubMed platform list information (http://www.ncbi.nlm. nih.gov/PubMed/) showed that out of 1.5 million works that deal with cancer analysis, only 5% are associated with mathematical modeling. However, it is clear that mathematics could contribute significantly to areas of experimental cancer analysis since there is currently a wealth of experimental information which needs a systematic analysis.

Even in these conditions, in the last decade, mathematical modeling and machine simulation of cancer has multiplied dramatically (e.g., reviews like [30-35]. A broad range of strategies were developed, specializing in one or additional aspects of cancer. For example, genetic instability, natural selection or interactions of individual cell with each other or the environment have been modeled using methods of cellular automata and agent-based modeling. These discrete methods have the disadvantage of being difficult to use when we deal with tumors of significant size. (see [36-38] for samples of cellular automata modeling and [39,40] for samples of agent-based modeling). In systems at larger scales, the neoplasic cell population is of the order of 10⁶ or more, making these discrete methods unfitted. For these situations, the continuum methods provide the best approach. Early work, as well as [41-43], used ODE to model cancer as a uniform population and partial differential equation models restricted to spherical geometries. To assess the stability of spherical tumors to asymmetric perturbations and to characterize the degree of aggression [31, 44-48] use linear and weakly nonlinear analysis. The interactions of a growth with the microenvironment, like stress-induced limitations to growth, are studied in [30,49-54]. For the sake of simplicity, most of the modeling has considered single-phase (e.g., single cell species) tumors. To provide an elaborate account of growth non-uniformity, [50, 55, 56] have been developed a mixture of models.

The results of morphology instabilities on each avascular and vascular solid neoplasm growth have been recently studied using non-linear modeling. With the help of boundary integral methods, Cristini et al [47] performed the first absolutely nonlinear simulations of a time model of neoplasm growth within the avascular and vascular growth stages with arbitrary boundaries. The model from [47] has been extended in 3D by Li et al. [48] via adaptive boundary integral technique. The inclusion of angiogenesis and extratumoral environment has been performed by Zheng et al., [57]. By developing and coupling a level set implementation with a hybrid continuum-discrete growing model originally developed by Anderson & Chaplain [58]they found that low-nutrient (e.g., hypoxic) conditions could lead to morphological instability. Their work served as a building block for recent studies of the impact of therapy on neoplasm growth [59] and for studies of morphological instability and invasion [60-62]. Macklin & Lowengrub used a ghost cell/level set technique for evolving interfaces to check neoplasm growth in heterogeneous tissue and additional studied neoplasm growth as a function of the microenvironment [63]. Wise et al. [64] and Frieboes et al. [65] have developed a diffuse interface implementation of solid neoplasm growth for the study of the evolution of multiple neoplasm cell species, that was used in [65] to model the 3-D vascularised growth of malignant gliomas (brain tumors).

In biological systems, the fractal structure of area in which cells act and differentiate is important for their organization and emergence of the hierarchical network of multiple cross-interacting cells, sensitive to external and internal conditions. The biological phenomena occur within the area whose dimensions aren't represented solely by integers (1, 2, 3, etc.) of Euclidean space. Particularly, malignant tumors [53-56] grow in a space with non-integer dimension, i.e. fractal dimension. The analytical formulae describing the time-dependence of the temporal fractal dimension and scaling reproduce the expansion of the Flexner–Jobling rat’s neoplasm in particular and growth of different rat’s tumors generally. The results of some calculations indicated that the formula derived for the time-dependent temporal fractal dimension and the scaling factor describe the experimental data obtained by Schrek for the Brown-Pearce rabbit’s neoplasm growth within the fractal time-space [3, 66-68].

In our assertion, fractal space(-time) consists in developing the consequences of the withdrawal of space(-time) differentiability’s hypothesis and acquiring a fractal geometry, namely space(-time) becomes explicitly dependent on the observation scale [69].

On the other hand, of great use in our further reasonings will be the fact that in many biological systems it is possible to empirically demonstrate the presence of attractors that operate starting from different initial conditions (Ivancevic). Some of these attractors are points, some are closed curves, while the others have non–integer, fractal dimension and are termed “strange attractors” [70]. It has been proposed that a prerequisite for proper simulating tumor growth by computer is to establish whether typical tumor growth patterns are fractal. The fractal dimension of tumor outlines was empirically determined using the box-counting method [71]. In particular, fractal analysis of a breast carcinoma was performed using a morphometric method, which is the box-counting method applied to the mammogram as well as to the histological section of a breast carcinoma [72].

If tumor growth is chaotic, this could explain the unreliability of treatment and prediction of tumor evolution. More importantly, if chaos is established, this could be used to adjust strategies for fighting cancer. Treatment could include some form of chaos control and/or anti-control.

4. A few words about holography

“Although it generates a three-dimensional image, a hologram is most often recorded on a photographic plate or a flat piece of film. Moreover, producing a hologram does not imply, in the conventional sense, the recording of an image. To better understand this apparent paradox and, as a result, the way holography works, we have to begin with the main principles.

In conventional imaging techniques, e.g. photography, what is being recorded is merely the intensity distribution in the original scene. Thus, all information about the optical paths to different parts of the scene is lost.

The unique property of holography is the method of recording both the phase and the amplitude of the light waves from an object. Since all recording materials respond only to the intensity in the image, it is mandatory to convert the phase information into intensity variations. Holography accomplishes this by using coherent illumination and introducing, as shown in Fig. 6, a reference beam derived from the same source. The photographic film records the interference pattern produced by this beam and the light waves scattered by the object in cause.

Since the intensity at any given point in this pattern of interference also depends on the phase of the object wave, the resulting recording (the hologram) contains information on the phase as well as the amplitude of the object wave. If the hologram is illuminated once again with the original reference wave, as shown in Fig. 7, it reconstructs the original object wave.

Figure 6.

Hologram recording: the interference pattern produced by the reference wave and the object wave is recorded.

An observer looking through the hologram sees a perfect three-dimensional image. This image exhibits all the effects of perspective, and depth of focus when photographed, that characterized the original object.

4.1. Early development

Gabor’s historical experiment of holographic imaging [73] consisted in a transparency formed of opaque lines on a clear background which was illuminated with a collimated beam of monochromatic light, the interference pattern produced by the directly transmitted beam (the reference wave) and the light scattered by the lines on the transparency being recorded on a photographic plate. When the hologram (i.e. a positive transparency made from this photographic negative) was illuminated with the original collimated beam, it produced two diffracted waves, one which reconstructed an image of the object in its original location, and the other, with identical amplitude but an opposite phase, which formed a second, conjugate image.

Figure 7.

Image reconstruction: light diffracted by the hologram reconstructs the object wave.

An important flaw of this method of image reconstruction was the poor quality of the resulting image, due to the fact that it was degraded by the conjugate image that was superimposing on it as well as by the scattered light from the directly transmitted beam.

Leith and Upatnieks [74-76] found a solution to the above-mentioned problem developing a off-axis reference beam technique presented schematically in Figs. 6 and 7. They used a separate reference wave incident on the photographic plate at an appreciable angle to the object wave. As a result, when the hologram was illuminated with the original reference beam, the two images were separated by large enough angles from the directly transmitted beam, and from each other, thus ensuring that the images not overlap.

The improvement of the off-axis technique, and, in equal measure, the invention of the laser, which provided a powerful source of coherent light, resulted in a surge of activity in holography that led to several crucial applications.

4.2. The in-line hologram

Let us now look upon the optical system presented in Fig. 8 in which the object (a transparency containing small opaque details on a clear background) is illuminated by a collimated beam of monochromatic light along an axis normal to the photographic plate.

Figure 8.

Optical system used to record an in-line hologram.

We can observe two components of the incident light. The first is the directly transmitted wave, which is a plane wave whose amplitude and phase do not vary across the photographic plate. Thus, its complex amplitude can be noted as a real constant r. The second one is a weak scattered wave whose complex amplitude at any point (x, y) on the photographic plate can be noted as o(x, y), where |o(x, y)| << r.

From these it can be shown that the resulting complex amplitude is the sum of these two complex amplitudes, and because of that the intensity at this point is

$$I(x,y)={\left|r+o(x,y)\right|}^2=r^2+{\left|o(x,y)\right|}^2+ro(x,y)+r{o}^{\ast }(x,y)$$E1

where o*(x, y) is the complex conjugate of o(x, y).

A ‘positive’ transparency (the hologram) is then made by contact printing from this recording. Therefore it can be assumed that this transparency is processed so that its amplitude transmittance (the ratio of the transmitted amplitude to that incident on it) can be written as

$$t=t_0+\beta TI$$E2

where t₀ is a constant background transmittance, T is the exposure time and β is a parameter determined by the photographic material used and the processing conditions, the amplitude transmittance of the hologram is

$$t(x,y)=t_0+\beta T[r^2+{\left|o(x,y)\right|}^2+ro(x,y)+r{o}^{\ast }(x,y)]$$E3

Then, the hologram is illuminated, as shown in Fig. 9, with the same collimated beam of monochromatic light employed to produce the original recording. Since the complex amplitude at any point in this beam is, aside from a constant factor, the same as that in the original reference beam, the complex amplitude transmitted by the hologram can be written as

$$u(x,y)=rt(x,y)=r(t_0+\beta T{r}^2)+\beta Tr{\left|o(x,y)\right|}^2+\beta T{r}^{2}o(x,y)+\beta T{r}^2{o}^{\ast }(x,y)$$E4

The right-hand side of (4) contains four terms. The first, r(t₀+βTr²), which represents a uniformly attenuated plane wave, corresponds to the directly transmitted beam.

The second, βTr |o(x, y)|², can be neglected, because is extremely small, compared to the other terms.

The third term, βTr²o(x, y), is, except for a constant factor, identical with the complex amplitude of the scattered wave from the object and has the property of reconstructing an image of the object in its original position. Due to the fact that this image forms behind the hologram, and the reconstructed wave appears to diverge from it, it is a virtual image.

The fourth term, βTr²o*(x, y), represents a wave similar to the object wave, but having an opposite curvature. This wave converges to form a real image (the conjugate image) at the same distance in front of the hologram.

With an in-line hologram, an observer viewing one image sees it superimposed on the out-of-focus twin image as well as a strong coherent background. Another problem is that the object must have a high average transmittance in order for the second term on the right-hand side of (4) to be negligible. Thus, it is possible to form images of fine opaque lines on a transparent background, but not vice versa. Finally, the hologram must be a ‘positive’ transparency. If the initial recording is used directly, β in (2) is negative, and the reconstructed image can be considered a photographic negative of the object.

Figure 9.

Optical system used to reconstruct the image with an in-line hologram, showing the formation of the twin images.

4.3. Off-axis holograms

In order to understand the formation of an image by an off-axis hologram, we must consider the recording arrangement shown in Fig. 10, in which (for simplicity) the reference beam is a collimated beam of uniform intensity, derived from the same source as the one used to illuminate the object.

Figure 10.

The off-axis hologram: `recording.

The complex amplitude at any point (x, y) on the photographic plate due to the reference beam can then be written as

$$r(x,y)=r\exp (i2\pi \xi x)$$E5

where ξ=(sinθ)/λ. Since only the phase of the reference beam varies across the photographic plat and because of the object beam, for which both the amplitude and phase vary, we can write the following:

$$o(x,y)=\left|o(x,y)\exp \left[-i\varphi (x,y)\right]\right|$$E6

The resultant intensity is, therefore,

$$\begin{array}{l}I(x,y)={\left|r(x,y)+o(x,y)\right|}^2={\left|r(x,y)\right|}^2+{\left|o(x,y)\right|}^2+\\ +r\left|o(x,y)\right|\exp \left[-i\varphi (x,y)\right]\exp (-i2\pi \xi x)+r\left|o(x,y)\right|\exp \left[i\varphi (x,y)\right]\exp (i2\pi \xi x)=\\ =r^2+{\left|o(x,y)\right|}^2+2r\left|o(x,y)\right|\cos [2\pi \xi x+\varphi (x,y)]\end{array}$$E7

The amplitude and phase of the object wave are encoded as amplitude and phase modulation, respectively, of a set of interference fringes equivalent to a carrier with a spatial frequency of ξ.

If, as in (2), we assume that the amplitude transmittance of the processed photographic plate is a linear function of the intensity, the resultant amplitude transmittance of the hologram is

$$\begin{array}{l}t(x,y)={t^{\prime }}_0+\beta T{\left|o(x,y)\right|}^2+\beta Tr\left|o(x,y)\right|\exp [-i\varphi (x,y)]\exp [-i2\pi \xi x]+\\ +\beta Tr\left|o(x,y)\right|\exp [i\varphi (x,y)]\exp [i2\pi \xi x]\end{array}$$E8

where $t_0^{\text{'}}$ = t₀ +βTr² is a constant background transmittance.

When the hologram is illuminated for a second time with the original reference beam, as shown in Fig. 11, the complex amplitude of the transmitted wave can be written as

$$\begin{array}{l}u(x,y)=r(x,y)t(x,y)={t^{\prime }}_{0}r\exp (i2\pi \xi x)+\beta Tr{\left|o(x,y)\right|}^2\exp (i2\pi \xi x)+\\ +\beta T{r}^{2}o(x,y)+\beta T{r}^2{o}^{\ast }(x,y)\exp (i4\pi \xi x)\end{array}$$E9

The first term on the right-hand side of (9) corresponds to the directly transmitted beam, while the second term generates a halo surrounding it, with approximately twice the angular spread of the object. The third term is identical to the original object wave, except for a constant factor βTr², and produces a virtual image of the object in its original position. The fourth term corresponds to the conjugate image which, in this case, is a real image. If the offset angle of the reference beam is taken large enough, the virtual image can be separated from the directly transmitted beam and the conjugate image.

In this setup, corresponding points on the real and virtual images are located at equal distances from the hologram, but on opposite sides of it. Because the depth of the real image is inverted, it is called a pseudoscopic image, as opposed to the normal, or orthoscopic, virtual image. It should also be mentioned that the phase of the reconstructed image is influenced only by the sign of β, always resulting in a “positive” image, even in the case of the hologram recording being a photographic negative.

Figure 11.

The off-axis hologram: image reconstruction.

4.4. Recording materials

Several recording materials have been used for holography [77]. Table 1 lists the principal characteristics of those that have been found most useful.

Material	Exposure; J/m^2	Resolution; mm^-1	Processing	Type	ηmax; (diffraction efficiency)
Photographic	≈1.5	≈ 5000	Normal	Amplitude	0.06
			Bleach	Phase	0.60
DCG; (dichromated gelatin)	10^2	10000	Wet	Phase	0.90
Photoresists	10^2	3000	Wet	Phase	0.30
Photopolymers	10-10^4	5000	Dry	Phase	0.90
PTP; (photothermoplastics)	10^-1	500-1200	Dry	Phase	0.30
BSO; ( Bi12SiO20 photorefractive crystals)	10	10000	None	Phase	0.20

Table 1.

Recording materials for holography

High-resolution photographic plates and films were the first materials used to record holograms. These are used widely even now, due to the fact that they exhibit relatively high sensitivity when compared to other hologram recording materials [78]. Moreover, they can be dye sensitized so that their spectral sensitivity matches the most commonly used laser wavelengths.

Combining the high sensitivity of photographic materials with the high diffraction efficiency, low scattering and high light-stability of DCG (dichromated gelatin) [79] was made possible by the silver-halide sensitized gelatin technique.

In positive photoresists, such as Shipley AZ-1350, the areas exposed to light become soluble and are washed away during development to produce a relief image [80].

Several organic materials can be activated by a photosensitizer to produce refractive index changes, because they suffer photopolymerization, when exposed to light [81]. A commercial photopolymer, coated on a polyester film base (DuPont OmniDex) that can be used to produce volume phase holograms with high diffraction efficiency is being currently produced [82].

Photothermoplastics (PTP) - a hologram can be recorded in a multilayer structure consisting of a glass or Mylar substrate coated with a thin, transparent, conducting layer of indium oxide, a photoconductor, and a thermoplastic [83,84].

When a photorefractive crystal is exposed to a spatially varying light pattern, electrons are liberated in the illuminated areas. These electrons migrate to adjacent dark regions, being trapped there. The spatially varying electric field produced by this space-charge pattern modulates the refractive index through the electro-optic effect, producing the equivalent of a phase grating. The space charge pattern can be removed by uniformly illuminating the crystal, after which another recording can take place [85,86].

It is essential to use coherent illumination for maximizing the visibility of the interference fringes formed by the object and reference beams, in the process of recording a hologram. In addition to being spatially coherent, the coherence length of the light must be much greater than the maximum value of the optical path difference between the object and the reference beams in the recording system. Lasers are, as a result, employed almost universally as light sources for recording holograms.” (The text in quotation marks was reproduced from Hariharan P. [165]).

Consequently, to get a hologram, one needs a laser, which provides a powerful source of coherent light, and a ‘recording material’ which records the interference pattern produced by a reference beam and the light waves scattered by the object.

5. Tumor-associated ECM or tumor microenvironment, nonlinear medium with holographic properties

5.1. Extracellular matrix and tumor microenvironment

Within tissue, cells are surrounded by a meshwork of proteins and proteoglycans collectively called the extracellular matrix (ECM), which compartmentalizes tissues. The ECM is divided into two distinct layers:

1. the basement membrane, which is composed of sheet-like layers of ECM and lies under epithelial cells segregating tissues into functionally distinct regions;

2. the interstitial matrix, which exists within intercellular space. The ECM serves multiple functions that are critical for embryonic development and wound repair. These functions include providing tissues with shape and flexibility and acting as a cushion to absorb external pressure. The ECM also serves as a base for cell anchorage, which mediates cell polarity, intercellular signaling, and assists in migration. The key to the ECM function lies in its unique composition and structure. The ECM is constructed in a specific pattern that is critical to its ability to carry out these functions and alterations in the expression level or arrangement of proteins within the ECM can be used to manipulate its function.

The most obvious function of the ECM is to provide structural support, shape, and stability for tissues. It does this by functioning as a base for cell anchorage. This base consists of three main structural components collagen, fibronectin, and elastic fibers, which bind to one another building a protein lattice upon which cells adhere.

Cell adherence to the ECM lattice provides cells support needed for cell migration. This is particularly important during embryonic development when cells are required to migrate into surrounding regions and differentiate into specific tissues [87]. A less obvious yet possibly more important function of the ECM in regards to tissue homeostasis and disease is its ability to mediate intracellular signaling. The ECM affects signaling through three main mechanisms:

1. cell – ECM interaction;

2. regulation of the bioavailability of growth factors;

3. the function of matricellular proteins. Cell attachment to the ECM via integrins induces signaling cascades that promote survival. Loss of cell-ECM contact can result in a form of apoptosis termed anoikis [88]. Anchorage-dependent survival is observed in most cells with the exception of red blood cells and inflammatory cells. However, tumor cells are often resistant to anoikis and can survive without a physical attachment to the ECM allowing them to successfully metastasize to distant tissues [89].

The ECM also affects cellular activity by serving as a reservoir for proteins required for proper tissue function and repair. This includes a plethora of growth factors and proteases. These pleiotropic molecules have been shown to robustly affect proliferation, survival and migration in numerous cell types. Once growth factors are secreted from cells, they often become embedded within the ECM and require ECM degradation by proteases such as elastase to release the active protein allowing it to interact with surrounding and transduce downstream signaling. The ability of the ECM to control the bioavailability of growth factors provides another means of regulating cellular activities and further explains how alterations in the makeup of the ECM as observed in diseases such as cancer affect cell response.

Matricellular proteins also reside in the ECM. They are a unique family of proteins that do not function as structural proteins but rather orchestrate the deposition of the ECM and mediate cell-cell and cell-ECM interactions. To do this, matricellular proteins interact directly with cell surface receptors, structural proteins, growth factors and proteases found within the ECM [90]. Their expression is found in every tissue begins early in development, persists throughout adulthood and is increased during tissue remodeling events. Matricellular proteins are critical regulators of many aspects of cell function including differentiation, survival, proliferation and migration making them necessary for proper tissue function. Not surprisingly, given their affect on cell-ECM mediated signaling pathways, matricellular proteins have been shown to strongly influence tumor growth.

For tumor cells to metastasize, the local ECM must be remodeled to create an environment conducive to tumor survival and progression. This includes altering the architecture and composition of the tumor-associated ECM or tumor microenvironment (TME) to facilitate tumor cell dissemination [91]. Changes in ECM architecture are primarily carried out by enzymes such as MMPs which assist in remodeling of the TME by degrading structural proteins such as collagen and fibronectin allowing tumor cells to freely navigate through the surrounding ECM. MMPs and other proteases assist in destruction of the first barrier tumor cells face to successful metastasis, the basement membrane. They degrade the underlying basement membrane allowing tumor cells to escape the primary tumor and invade into surrounding non-neoplasic tissues. MMPs continue to breakdown barriers in the surrounding ECM clearing a path to blood vessels where tumor cells will intravasate into the circulatory system and seed secondary tumors [92]. Destruction of the ECM by proteases also promotes tumor progression by facilitating the release of angiogenic and mitogenic factors bound within the ECM [93]. In a surprising unexpected twist, studies revealed that the breakdown of ECM proteins by MMPs was more complex than anticipated. In fact it is a highly organized process which results in the generation of both protumor and antitumor cleavage products [94].

Presence of the ECM is required for cellular survival therefore increased degradation of the ECM within the TME must be balanced by an increase in ECM synthesis. The development of a tumor, much like a wound, provokes a robust inflammatory response causing an influx of mast cells, macrophages and neutrophiles into the TME [95].

We may summarize that the extracellular matrix generates signaling cues that regulate cell behavior and orchestrate functions of cells in tissue formation and homeostasis. Microenvironmental signaling, a process that determines cell shape, motility, growth, survival and differentiation is highly influenced by the ECM properties: composition, three-dimensional organization and proteolytic remodeling. Recent studies have shown that misregulation of cell–ECM interactions can contribute to many diseases, including developmental, immune, haemostasis, degenerative and malignant disorders.

Consequently, the structure and the behavior of the tumor-associated ECM allows us to think of it as a non-differential medium, and as will be shown below, a medium which holds the properties of a hologram (capacity to memorize, interference abilities) and may become a source of forces. In other words, ECM and TME are very suitable candidates for a ‘recording material’.

5.2. Tumor-associated ECM as a non-differential fractal medium

We can simplify the dynamics of a biological system supposing that the motions on ECM take place on continuous but non-differentiable curves, i.e. fractal curves (for example, the Peano curve, the Koch curve or the Weierstrass curve [69,96,97].

Once this hypothesis has been accepted, some consequences of non-differentiability by SRT are evident [69,96]: i) the physical quantities that are used in describing the biological system dynamics are fractal functions, i.e. functions dependent both on spatial coordinates and time as well as on the scale resolution, δt/τ (identified here with dt/τ by means of the substitution principle [69,96]. We mention that in the standard biophysics, the physical quantities describing the dynamics of a biological system are continuous, but differentiable functions depending only on spatial coordinates and time; ii) the dynamics of the biological systems are given by the fractal operator $\widehat{d}$ /dt [98]:

$$\frac{\widehat{d}}{dt}=\frac{\partial }{\partial t}+\widehat{V}\cdot \nabla -i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\Delta$$E10

where

$$\widehat{V}=V_D-i{V}_F$$E11

is the complex velocity, V_D is the differentiable and resolution scale independent velocity, V_F is the non-differentiable and resolution scale dependent velocity, $\widehat{V}$ ∇ is the convective term,

$$\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\Delta =\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}\right)$$E12

is the dissipative term, D_F is the fractal dimension of the movement curve, λ is the space scale, τ is the time scale and λ²/τ is a coefficient specific to the fractal – non - fractal transition. For D_F any definition can be used (the Hausdorff – Besikovici fractal dimension, the Kolmogorov fractal dimension, etc. [97], but once such definition is accepted for D_F, it has to remain constant over the entire analysis of the complex fluid dynamics. In a particular case, for motions on Peano curves, D_F = 2 [97] of the complex fluid entities, the fractal operator (1) is reduced to Nottale’s operator ($\widehat{d}$ /dt)_Ν

$$\frac{\widehat{d}}{dt}=\frac{\partial }{\partial t}+\widehat{V}\cdot \nabla -i{D}_N\Delta$$

where D_N = λ²/τ is the Nottale’s coefficient associated to the fractal-non-fractal transition.

Applying the fractal operator (10) to the complex velocity (11) and accepting the principle of scale covariance [69,96] in the form:

$$\frac{\widehat{d}\widehat{V}}{dt}=-\nabla U$$E13

we obtain the motion equation:

$$\frac{\widehat{d}\widehat{V}}{dt}=\frac{\partial \widehat{V}}{\partial t}+\left(\widehat{V}\cdot \nabla \right)\widehat{V}-i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\Delta \widehat{V}=-\nabla U$$E14

where U is an external scalar potential. Equation (14) is a Navier – Stokes type equation. It means that at any point of a fractal path, the local acceleration term, ∂_t $\widehat{V}$_, the non-linearly (convective) term, ($\widehat{V}$ ∇) $\widehat{V}$, the dissipative term, $$\left({\lambda }^2/\tau \right){\left(dt/\tau \right)}^{\left(\frac{2}{D_F}\right)-1}\Delta \widehat{V}$$, and the external free term ∇U make their balance. Therefore, the biological fluid is assimilated to a “rheological” fractal fluid, whose dynamics are described by the complex velocities field, $\widehat{V}$, and by the imaginary viscosity type coefficient, $$i\left({\lambda }^2/\tau \right){\left(dt/\tau \right)}^{\left(\frac{2}{D_F}\right)-1}$$. The “rheology” of the fractal fluid can provide hysteretic properties to the biological fluid (the fractal fluid has a hysteresis cycle, memory, etc. [98-100].

For irrotational motions of the biological system entities

$$\begin{array}{l}\nabla \times \widehat{V}=0\text{,}\text{a}\\ \nabla \times {\widehat{V}}_D=0,\text{b}\\ \nabla \times {\widehat{V}}_F=0\text{c}\end{array}$$E15

we can choose $$\widehat{V}$$ of the form

$$\widehat{V}=-i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\nabla \ln \psi$$E16

where ϕ ≡ lnψ is the velocity scalar potential. By substituting (16) in (14) and using the method described in [98-100], it results

$$\frac{\widehat{d}\widehat{V}}{dt}=-i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\nabla \left[\frac{\partial \ln \psi }{\partial t}-i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\frac{\nabla \psi }{\psi }+U\right]=0$$E17

This equation can be integrated in a universal way and yields

$$\frac{{\lambda }^4}{{\tau }^2}{\left(\frac{dt}{\tau }\right)}^{\left(\frac{4}{D_F}\right)-2}\Delta \psi +i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\frac{\partial \psi }{\partial t}-\frac{U}{2}\psi =0$$E18

up to an arbitrary phase factor which may be set to zero by an appropriate selection of the phase of ψ. Relation (18) is a Schrödinger type equation. For motions on Peano curves, D_F = 2 [97] at Compton scale, which implies λ² /τ = ħ /2m₀ [69,96], with ħ the reduced Plank constant and m₀ the rest mass of the biological entities, the relation (18) becomes the standard Schrödinger equation:

$$\frac{{\hslash }^2}{2m_0}\Delta \psi +i\hslash \frac{\partial \psi }{\partial t}-\frac{U}{2m_0}\psi =0$$

If ψ = √ρ e^iS, with √ρ the amplitude and S the phase of ψ, the complex velocity field (16) takes the explicit form:

$$\begin{array}{l}\widehat{V}=\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\nabla S-i\frac{{\lambda }^2}{2\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\nabla \ln \rho \text{a}\\ V_D=\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\nabla S\text{b}\\ V_F=\frac{{\lambda }^2}{2\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\nabla \ln \rho \text{c}\end{array}$$E19

By substituting (19a-c) in (14) and separating the real and the imaginary parts, up to an arbitrary phase factor which may be set to zero by appropriate selection of the phase of ψ, we obtain:

$$\frac{\partial V_D}{\partial t}+\left(V_D\cdot \nabla \right)V_D=-\nabla \left(Q+U\right)$$E20

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho V_D\right)=0$$E21

with Q the specific fractal potential

$$Q=-2\frac{{\lambda }^4}{{\tau }^2}{\left(\frac{dt}{\tau }\right)}^{\left(\frac{4}{D_F}\right)-2}\frac{\Delta \sqrt{\rho }}{\sqrt{\rho }}=-\frac{V_F^2}{2}-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(\frac{2}{D_F}\right)-1}\nabla \cdot V_F$$E22

Equation (20) represents the specific momentum conservation law, while equation (21) represents the states density conservation law. By means of the fractal velocity, V_F, the specific fractal potential Q is a measure of non-differentiability of the biological entities trajectories, i.e. of their chaoticity. The equations (20)-(22) define the fractal hydrodynamics model (FHM). In such a context, the biological system can be considered a fractal fluid.

Thus, it can be concluded that: i)Any entity of the biological system is in a permanent interaction with the fractal medium by means of the specific fractal potential; ii) The fractal medium is identified with a non-relativistic fractal fluid described by equations (20)-(22); iii) For motions on Peano curves at Compton scale [69,96,97], the FHM reduces to a quantum hydrodynamic model (QHM). Indeed, according to our previous considerations the relations (19a-c) become

$$\begin{array}{l}\widehat{V}=\frac{\hslash }{m_0}\nabla S-i\frac{\hslash }{2m_0}\nabla \ln \rho \\ V_D=\frac{\hslash }{m_0}\nabla S\\ V_F=\frac{\hslash }{m_0}\nabla \ln \rho \end{array}$$

in order that the momentum and density conservation laws are given by (20) and (21), respectively, with V_D and V_F previously defined, and the specific fractal potential by the expression

$$Q=-\frac{{\hslash }^2}{2m_0}\frac{{\nabla }^2\sqrt{\rho }}{\sqrt{\rho }}=-\frac{V_F^2}{2}-\frac{\hslash }{2}\nabla \cdot V_F$$

Moreover the fractal medium is assimilated to Bohm subquantum level [96]; iv) The fractal velocity V_F cannot be regarded as actual mechanical motion; it contributes to the transfer of the specific momentum and the concentration of energy. This fact can easily be deduced from the absence of V_F in the states density conservation law, and from its role in the variational principle. Any interpretation of Q should take into account the “self” or internal nature of the specific momentum transfer. While the energy is stored in the form of mass motion and potential energy (as it is classically), a part of it is available elsewhere and only the total is conserved. Reversibility and the existence of eigenstates is ensured by the conservation of energy and specific momentum, but this also means that a Brownian motion [97] form of interaction with an external medium is denied; v) For Peano curves motions [96,97], at spatial scales higher than the dimension of the boundary layer and at temporal scales higher than the oscillation periods of the pulsating velocities which overlaps the average velocity of the biological fluid motions (for details see [101-103], the FHM reduces to the standard hydrodynamical model [104]; vi) Since the position vector of the biological system entity is assimilated with a stochastic Wiener type process [96,97], ψ is not only the scalar potential of a complex velocity (through ϕ ≡ lnψ) in the frame of FHM, but represents also the states density (through ψ²) in the frame of a Schrödinger type model. Thus it can be seen that the formalism of the FHM and the one of Schrödinger type are equivalent. In addition, the chaoticity, either by means of turbulence in the fractal hydrodynamics approach, or by means of stochasticization in the Schrödinger type approach, is generated only by the non-differentiability of the movement trajectories in a fractal space; vii) In the standard model (Landau’s scenario [104]) the Fourier spectrum is always discrete and cannot approximate a continuum spectrum that in case of a large number of frequencies will generate an unlimited number of spectral components as a result of their beats which appear due to the presence of nonlinearities in the biological fluid. Still, taking into account the standard model, the flow can never be exactly chaotic because, in case of multiple periodic functions, correlations tend to be null, although having an oscillating character. As a result, the transition towards chaotic behavior can be described by Landau’s scenario only in a biological system with an infinite number of degrees of freedom. In our case, when δ t/τ $\to$ 0 for D_F ≠ 2 the physical quantities that describe the dynamics of the biological system are no longer defined. So, in this approximation, a simulation of a system with an infinite number of degrees of freedom is used. Moreover, the possibility of the dynamic states generation should be noted, which is characterized by windows of regular oscillations interrupted by chaotic bursts, the transition between the two states being spontaneous, unpredictable and independent of any of the control parameters variation (turbulence through intermittency); viii) The fractal medium and in particular the subquantum level has some computational properties: viii1) bistability, which implies the existence of its fractality and in particular, for motions on Peano curves at Compton scales, of the quantum bit. And from here, the entire fractal logics and in particular the quantum one; viii2) the self-replication, which implies the existence of some specific self-copying mechanisms; viii3) memory, which implies hysteresis type mechanisms; viii4) self-similarity, which implies the holographic type behavior; viii5) polarization, which implies mechanisms of changing the “computational state” from a given to a desired one; viii6) depositing and transmitting the information etc. For details see [105].

6. Tumor self-seeding by CTC and hypoxia support the idea of complete holography

6.1. The self-seeding hypothesis of tumor growth

The unfolding of cancer cells from their original sites to different ones within the body, i.e metastasis, has, for many years, been regarded as a unidirectional journey. However some researchers conjointly consider that metastatic cancer cells can increase primary tumor growth, this fact being crucial for the planning and type of the cancer treatment.

The concept of growing self-metastasis, or tumor “self-seeding,” was first introduced at Memorial Sloan-Kettering Cancer Center, in a range of studies conducted by Drs. Joan Massagué, head of the Metastasis research facility, and Larry Norton, deputy physician-in-chief of the center’s breast cancer programs. In the studies conducted on mice, Dr. Massagué discovered that breast tumors express genes related to metastasis were growing quicker than tumors that didn't express these genes, even if the genes had no apparent role in increased cellular division or decreased cell death (Fig. 12). These results did not fit within the standard tumor growth theories. In 2006, the researchers theorized that cells that become independent from a tumor and colonize distant tissues may also return home to the microenvironment within which they initial developed via the cardiovascular system [106]. They tested their hypothesis in a mouse model of cancer and revealed their findings in 2009 in Cell [107].

In one particular experiment, they selected a non-metastatic breast cancer cell line and an isolated set of daughter cells from that line that had gained the ability over time to metastasize to the lungs. Consequently, they implanted the parent cells in one mammary gland and the metastatic daughter cells in the opposite gland to serve as “donor tumors”. They noticed that the daughter cells migrated to the lungs and to the tumors that were being formed by the parent cells in the opposite organ, accounting for 5 to 30% of the size of the parent tumors. Also, it was obvious that parent tumors seeded by daughter cells grew quicker than parent tumors that were implanted without daughter cells within the opposite gland.

This specific seeding behavior with daughter cells that spread to the bones and brain was noticed in the studies of colon cancer and skin cancer cell lines, but not when non-metastatic daughter cells were transplanted.

Furthermore, in various related follow-up laboratory experiments, the researchers demonstrated that cells from primary tumors can attract circulating metastatic tumor cells, and found several proteins that probably encourage this migration. They also found that the “come-back” metastatic cells mainly influenced primary tumor growth through the release of proteins that modify the tumor microenvironment, as well as blood vessels and immune cells.

Their hypothesis started to obtain support from different researchers at the United Nations agency. In 2009, Dr. Philip Hahnfeldt and his colleagues published the results of computer modeling studies conceived to search at the intersection of two biological phenomena found in tumors [108]. One of these phenomena consists in a small population of cancer cells acting like stem cells; thus, they could possess the ability to reproduce an infinite number of times, and also to generate secondary cancer cells that, with time, lose the ability to divide. The second phenomena is that tumor growth is restricted by the available space for growing. Normally, healthy cells are spatially separated by a space that is not available in the tumor, due to the fact that cancer cells grow tightly in a dense mass till all the available has been occupied and the cellular division stops. At the periphery of the tumor, where the normal tissues density decreases, cancer cells continue to multiply and expand, increasing the size of the tumor.

Figure 12.

In the self-seeding concept of cancer growth and metastasis, a mobile tumor cell can take one of five different pathways in the body. A – evade and return to the primary tumor, using only the close ECM and not the systemic circulation; B – escape into the systemic circulation and then return to the original tumor; C – migrate through the systemic circulation and grow a metastatic tumor elsewhere in the organism; D – evade and return to the metastatic tumor, not using the systemic circulation; E –escape and return to the metastatic tumor through the systemic circulation.

Their models showed a crucial relation between cell migration, cell death, and tumor growth. When the offspring of a cancer stem cell within the model did not migrate or die spontaneously, tumor growth remained constant at around 110 cells. On the other hand, a high death rate among the non-stem cell progeny combined with a high cell migration rate produced the biggest tumors within the shortest amount of your time, to virtually 100,000 cells in over three years.

This theoretical phenomena − accelerated growth jump-started by a high rate of growth death − has potential implications for the clinical treatment of cancer. Traditional cytotoxic therapy medicine kill massive numbers of speedily dividing cancer cells, however might not have an effect on cancer stem cells in each tumor type.

In the light of the above, we think of the CTC returning to the initial tumor site and fueling the primary tumor growth or even grow a new tumor as a particular case of complete holography (i.e. a hologram which does not represent only the virtual object’s image, but becomes the very object - which we believe, is a characteristic of the living organisms).

6.2. Hypoxia and cancer

Vascularized tissues is the trigger factor for a large number of cellular processes, combined with an adaptive response. [109,110]. Following the drop in oxygen supplies, cells start to adapt to the less favorable environment and to initiate a vascularization| process in order for them to raise the local oxygen supply. In the center of the hypoxic response is the angiogenic shift, with production of potent angiogenic factors such as vascular endothelial growth factor (VEGF). Hypoxia, which is present in many solid malignancies, means that oxygen level in tumors corresponds to around 1.5% [109]. This is caused on one hand by the result of the abnormal vascularisation in tumors, that is short in supply oxygen to the sometimes rapidly expanding malignant lesion and on the other hand, the existence of areas with acute lack of oxygen, resulting in necroses and cell death on a large scale. In breast cancer these forms are usually related to clinically aggressive behavior. Markers for hypoxia like HIF-1a have been connected to extremely malignant features and could be relevant prognostic markers for distinguishing subgroups of breast cancer with certain malignant properties [111,112]. There is a current discussion whether or not hypoxia contributes to increase the aggressiveness of tumors or if aggressive tumors have more widespread hypoxia, but, apparently, one explanation doesn't essentially exclude the opposite. Recent analysis has shown that ‘the hypoxia response’ in tumors may be used to conceive new treatment methods [113]. Emerging cancer therapies will most definitely put more focus on specific targeting of hypoxic processes.

Human cancers are characterized by intratumoral hypoxia that results from the proliferation of deregulated cell and the physiological responses that is triggered by it have impact on all aspects of cancer progression, together with im mortalization, transformation, differentiation, genetic instability, ontogeny, metabolic adaptation, autocrine protein communication, invasion, metastasis, and resistance to therapy.

We assume the relationship between hypoxia and aggressive tumors may be due to the presence of the coherent wave laser with oxygen of metastatic tumor cells in the area, where the produced oxygen gradients lead to oxygen consumption. It has been already shown that laser photocoagulation is effective in the treatment of diabetic retinopathy, in a series of major studies [114-117]. The oxygen-consumption may be based on a multilayer solution to Fick’s law of diffusion, yet the essence is that the oxygen consumption is greatest where the oxygen gradient changes most rapidly [118-120].

All the above considerations and hypoxia’s impact on all critical aspects of cancer progression support the idea that, the metastatic tumor cells moving through the systemic circulation (and not necessarily in there), may be considered a travelling wave chemically pumped type laser with oxygen.

7. Basic model

7.1. The PDE cancer-invasion model

We consider and present in what follows in extenso, the basic mathematical model of growth of a generic solid tumor, which is assumed just been vascularised, i.e. a blood supply has been established. Let us focus on four key variables involved in tumor cell invasion, in order to produce a minimal model, namely tumor cell density (denoted by n), matrix-degradative enzymes (MDE) concentration (denoted by m), the complex mixture of macromolecules from the extracellular material’s (MM) concentration (denoted by f) and the oxygen concentration (denoted by c). Each of the four variables (n, m, f, c) is a function of the spatial variable x and time t. Firstly, we have to define a system of coupled non-linear partial differential equations to model tumor invasion of surrounding tissue.

We make the assumption that the ECM is a mixture of MM (e.g. collagen, fibronectin, laminin and vitronectin) only and not any other cells. Most of the MM of the ECM which are important for cell adhesion, spreading and motility are fixed or bound to the surrounding tissue. MDEs are important at many stages of tumor growth, invasion and metastasis, and they interact with inhibitors, growth factors and tumor cells in a very complex way. Yet it is widely accepted that tumor cells produce MDEs which locally degrade the ECM. As well as creating space into which tumor cells may be transported by simple diffusion (random motility), we can assume that this also results in a gradient of these bound cell-adhesion molecules, such as fibronectin. As a result, while the ECM may be a barrier to normal cell movement, it also represents a substrate to which cells may adhere and move upon. The presence of a minimum of ECM elements is a requirement for the growth and survival of most mammalian cells, and indeed these cell will migrate up a gradient of bound (i.e. non-diffusible) cell-adhesion molecules in the in vitro cultures [121-126].

We can define haptotaxis as a directed migratory response of the cells to gradients of fixed or bound chemicals (i.e. non-diffusible chemicals). While studies have not yet clearly shown haptotaxis occurs in an in vivo situation, given the structure of human tissue, it is not without reason to assume that haptotaxis is a major component of directed movement in tumor cell invasion. Indeed, there has been much recent effort to characterize such directed movement [125-127]. We therefore will treat this directed movement of tumor cells in this model as haptotaxis, i.e. a response to gradients of bound MM such as fibronectin. To incorporate this response in the mathematical model, we take the haptotactic flux to be J_hapto = χ n ∇ f, where χ > 0 is the (constant) haptotactic coefficient.

As we stated early, the only other contribution to tumor cell motility in this model is the random motion. To describe the random motility of the tumor cells, we assume a flux of the form J_rand = −D_n∇n, where D_n is the constant random motility coefficient.

We only model the tumor cell migration at this level, as all other tumour cell processes, such as proliferation, adhesion and death will be treated at a single cell level within the hybrid discrete-continuum model. The conservation equation for the tumour cell density n can therefore be written as

$$\frac{\partial n}{\partial t}+\nabla \left(J_{rand}+J_{hapto}\right)=0$$

and hence the partial differential equation governing tumor cell motion (in the absence of cell proliferation) is

$$\frac{\partial n}{\partial t}=D_n{\nabla }^{2}n-\chi \nabla \left(n\nabla f\right)$$E23

The ECM is known to contain many MM, including fibronectin, laminin and collagen, which can be degraded by MDEs [128,129]. We assume that the MDEs degrade ECM upon contact and hence the degradation process is modeled by the following simple equation

$$\frac{\partial f}{\partial t}=-\delta mf$$E24

where δ is a positive constant.

Active MDEs are produced (or activated) by the tumor cells, diffuse throughout the tissue and undergo some form of decay (either passive or active). The equation governing the evolution of MDE concentration is therefore given by

$$\frac{\partial m}{\partial t}=D_m{\nabla }^{2}m+g(n,m)-h(n,m,f)$$E25

where D_m is a positive constant, the MDE diffusion coefficient, g is a function modeling the production of active MDEs by the tumor cells and h is a function modeling the MDE decay. For simplicity we assume that there is a linear relationship between the density of tumor cells and the level of active MDEs in the surrounding tissues (not taking into consideration the amount of enzyme precursors secreted and the presence of endogenous inhibitors) and so these functions will be g = μ n (MDE production by the tumor cells) and h = λ m (natural decay), respectively.

The fact solid tumors need oxygen to grow and invade is a well-known one. Oxygen is assumed to diffuse into the MM, decay naturally and be consumed by the tumor. We assume that oxygen production is proportional to the MM density. This is a crude way of modeling an angiogenic oxygen supply for a more appropriate way of modeling the angiogenic network. The oxygen equation then has the form,

$$\frac{\partial c}{\partial t}=D_c{\nabla }^{2}c+\beta f-\gamma n-\alpha c$$E26

where D_c, β, γ, α are positive constants representing the oxygen diffusion coefficient, production, uptake and natural decay rates, respectively.

The complete system of equations describing the interactions of the tumor cells, MM, MDEs and oxygen as detailed above, is

$$\begin{array}{l}\frac{\partial n}{\partial t}=\stackrel{\text{random motility}}{\overbrace{D_n{\nabla }^{2}n}}-\stackrel{\text{haptotaxis}}{\overbrace{\chi \nabla \left(n\nabla f\right)}}\text{a}\\ \frac{\partial f}{\partial t}=-\stackrel{\text{degradation}}{\overbrace{\delta mf}}\text{b}\\ \frac{\partial m}{\partial t}=\stackrel{\text{diffusion}}{\overbrace{D_m{\nabla }^{2}m}}+\stackrel{\text{production}}{\overbrace{\mu n}}-\stackrel{\text{decay}}{\overbrace{\lambda m}}\text{c}\\ \frac{\partial c}{\partial t}=\stackrel{\text{diffusion}}{\overbrace{D_c{\nabla }^{2}c}}+\stackrel{\text{production}}{\overbrace{\beta f}}-\stackrel{\text{uptake}}{\overbrace{\gamma n}}-\stackrel{\text{decay}}{\overbrace{\alpha c}}\text{d}\end{array}$$E27

where D_n, D_m and D_c are the tumor cell, MDE and oxygen diffusion coefficients, respectively, χ is the haptotaxis coefficient and δ, μ, λ, β, γ and α are positive constants. We should also note that cell–matrix adhesion is modeled here by the use of haptotaxis in the cell equation, i.e. directed movement up gradients of MM. Therefore, χ maybe considered as relating to the strength of the cell–matrix adhesion.

7.2. The PDE cancer-invasion model via scale relativity theory

The presence of the fractal medium implies the substitution of the standard derivative d/dt with the fractal operator (10). Then the system (27a-d) becomes

$$\begin{array}{l}\frac{\widehat{d}n}{dt}=\frac{\partial n}{\partial t}+\left(\widehat{V}\cdot \nabla \right)n-i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\Delta n=D_n{\nabla }^{2}n-\chi \nabla \left(n\nabla f\right)\text{a}\\ \frac{\widehat{d}f}{dt}=\frac{\partial f}{\partial t}+\left(\widehat{V}\cdot \nabla \right)f-i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\Delta f=-\delta mf\text{b}\\ \frac{\widehat{d}m}{dt}=\frac{\partial m}{\partial t}+\left(\widehat{V}\cdot \nabla \right)m-i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\Delta m=D_m{\nabla }^{2}m+\mu n-\lambda m\text{c}\\ \frac{\widehat{d}c}{dt}=\frac{\partial c}{\partial t}+\left(\widehat{V}\cdot \nabla \right)c-i\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\Delta c=D_c{\nabla }^{2}c+\beta f-\gamma n-\alpha c\text{d}\end{array}$$E28

or more explicitly, by separating the scales of interaction, for the differentiable scale

$$\begin{array}{l}\frac{\partial n}{\partial t}+\left({\widehat{V}}_D\cdot \nabla \right)n=D_n{\nabla }^{2}n-\chi \nabla \left(n\nabla f\right)\text{a}\\ \frac{\partial f}{\partial t}+\left({\widehat{V}}_D\cdot \nabla \right)f=-\delta mf\text{b}\\ \frac{\partial m}{\partial t}+\left({\widehat{V}}_D\cdot \nabla \right)m=D_m{\nabla }^{2}m+\mu n-\lambda m\text{c}\\ \frac{\partial c}{\partial t}+\left({\widehat{V}}_D\cdot \nabla \right)c=D_c{\nabla }^{2}c+\beta f-\gamma n-\alpha c\text{d}\end{array}$$E29

and for the fractal scale

$$\begin{array}{l}\left({\widehat{V}}_F\cdot \nabla \right)n=\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\Delta n\text{a}\\ \left({\widehat{V}}_F\cdot \nabla \right)f=\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\Delta f\text{b}\\ \left({\widehat{V}}_F\cdot \nabla \right)m=\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\Delta m\text{c}\\ \left({\widehat{V}}_F\cdot \nabla \right)c=\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\Delta c\text{d}\end{array}$$E30

Thus, the transport equations (27a-d) are generalized by involving the convective terms $$\left({\widehat{V}}_F\cdot \nabla \right)n,\left({\widehat{V}}_F\cdot \nabla \right)f,\left({\widehat{V}}_F\cdot \nabla \right)m,\left({\widehat{V}}_F\cdot \nabla \right)c$$ at differentiable scale. Moreover, at the fractal level one specifies new transport mechanisms where the convective effects are balanced by dissipative ones.

Now, the transport equations for the fractal to non-fractal transition are obtained by substracting the relations (29a) and (30a), (29b) and (30b), (29c) and (30c), (29d) and (30d) and using the substitution V=V_D-V_F. One gets

$$\begin{array}{l}\frac{\partial n}{\partial t}+\left(\widehat{V}\cdot \nabla \right)n=\left[D_n-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\right]\Delta n-\chi \nabla \left(n\nabla f\right)\text{a}\\ \frac{\partial f}{\partial t}+\left(\widehat{V}\cdot \nabla \right)f=-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\Delta f-\delta mf\text{b}\\ \frac{\partial m}{\partial t}+\left(\widehat{V}\cdot \nabla \right)m=\left[D_m-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\right]\Delta m+\mu n-\lambda m\text{c}\\ \frac{\partial c}{\partial t}+\left(\widehat{V}\cdot \nabla \right)c=\left[D_c-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\right]\Delta c+\beta f-\gamma n-\alpha c\text{d}\end{array}$$E31

Assuming now both the coherence fractal to non-fractal and harmonic type behavior for the f field the system of equations (31a-d) becomes

$$\begin{array}{l}\frac{\partial n}{\partial t}=\left[D_n-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\right]\Delta n-\chi \nabla \left(n\nabla f\right)\text{a}\\ \frac{\partial f}{\partial t}=-\delta mf\text{b}\\ \frac{\partial m}{\partial t}=\left[D_m-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\right]\Delta m+\mu n-\lambda m\text{c}\\ \frac{\partial c}{\partial t}=\left[D_c-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}\right]\Delta c+\beta f-\gamma n-\alpha c\text{d}\end{array}$$E32

7.3. Non-dimensionalization and parameterization

For us to utilize realistic parameter values, we must first non-dimensionalize the equations in the standard formalism. We therefore rescale the distance with an appropriate length scale L (e.g. the maximum invasion distance of the cancer cells at the first stage of invasion, approximately 1 cm), time with τ (e.g. the average time taken for mitosis to occur, approximately 8...24 h [130], tumor cell density with n₀, ECM density with f₀, MDE concentration with m₀ and oxygen concentration with c₀ (where n₀, f₀, m₀ and c₀ are appropriate reference variables). Therefore, setting

$$\tilde{n}=\frac{n}{n_0},\text{ <mover accent="true">f <mo>˜</mo></mover>}=\frac{f}{f_0},\text{ <mover accent="true">m <mo>˜</mo></mover>}=\frac{m}{m_0},\text{ <mover accent="true">c <mo>˜</mo></mover>}=\frac{c}{c_0},\tilde{x}=\frac{x}{L},\tilde{t}=\frac{t}{t_0}$$

in (27) and dropping the tildes for notational convenience, we obtain the scaled system of equations

$$\begin{array}{l}\frac{\partial n}{\partial t}=\stackrel{\text{random motility}}{\overbrace{d_n{\nabla }^{2}n}}-\stackrel{\text{haptotaxis}}{\overbrace{\rho \nabla \left(n\nabla f\right)}}\text{a}\\ \frac{\partial f}{\partial t}=-\stackrel{\text{degradation}}{\overbrace{\eta mf}}\text{b}\\ \frac{\partial m}{\partial t}=\stackrel{\text{diffusion}}{\overbrace{d_m{\nabla }^{2}m}}+\stackrel{\text{production}}{\overbrace{kn}}-\stackrel{\text{decay}}{\overbrace{\sigma m}}\text{c}\\ \frac{\partial c}{\partial t}=\stackrel{\text{diffusion}}{\overbrace{d_c{\nabla }^{2}c}}+\stackrel{\text{production}}{\overbrace{\nu f}}-\stackrel{\text{uptake}}{\overbrace{\omega n}}-\stackrel{\text{decay}}{\overbrace{\phi c}}\text{d}\end{array}$$E33

where d_n= τ $$D_n-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}$$ /L², ρ=τχ f₀/L², η=τm₀δ, d_m=τ $$D_m-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}$$ /L², κ=τμn₀/m₀, σ=τλ, d_c =τ $$D_c-\frac{{\lambda }^2}{\tau }{\left(\frac{dt}{\tau }\right)}^{\left(2/D_F\right)-1}$$ /L², ν=τ f₀β/c₀, ω=τn₀γ/c₀, ϕ=τα.

The cell cycle time can be highly variable (particularly the G1 phase) and in fact depends on the specific tumor taken under consideration. As an approximate reference time we take τ = 16 h, halfway between 8...24 h [130]. The cell motility parameter D_n ~ 10⁻⁹ cm² s⁻¹ was estimated from available experimental evidence [131]. Tumor cell diameters again will vary depending on the type of tumor being considered but are in the range 10...100 μm [132] with an approximate volume of 10⁻⁹ to 3 × 10⁻⁸ cm³ [133,134]. We will assume that a tumor cell has the volume 1.5 × 10⁻⁸ cm³ and therefore take n₀ = 6.7 × 10⁷ cells cm⁻³. The haptotactic parameter χ ~ 2600 cm² s⁻¹ M⁻¹ was estimated to be in line with that calculated in [135] and the parameter f₀ ~ 10⁻⁸ to 10⁻¹¹ M was taken from the experiments of [136]. We took D_m to be 10⁻⁹ cm² s⁻¹, which is perchance small for a diffusing chemical, but current studies imply that it is in fact a combination of the MDE and MM, and, as a result, the MM degrades and diffuses very little [137]. An in vivo estimate for the MDE concentration m₀ is rather difficult to obtain since no value (that we are aware of) has been currently determined and we also know that certain inhibitors (e.g. tissue inhibiting metalloproteases) are produced within the ECM which affects the MDE concentration. Plasma levels of specific MDEs have been measured (e.g. MMP-2 [138]) and are approximately 130 ng ml⁻¹ with further increases observed in patients with cancer [139]. How is this related to the MDE concentration within the ECM is not clear and we have therefore left this parameter undefined. Estimates for the kinetic parameters μ, λ and δ were not available since these are rather hard to obtain experimentally – and thus we use the values of [135]. The diffusion rate of oxygen through water is D_c = 10⁻⁵ cm² s⁻¹ and also, the oxygen consume rate of the cells is 6.25 × 10⁻¹⁷ M cells⁻¹ s⁻¹ [134]. The background oxygen concentration estimation within the tissue was somehow difficult to be done as it depends on the tissue vascularization. If we set the value of the oxygen concentration in the blood supplying the tumor/tissue to be 0.15 ml O₂ per ml of blood, since we know that 1 M of oxygen occupies 22400 ml then there is 0.15/22400 M O₂ ml⁻¹ = 6.7 x 10⁻⁶ M O₂ ml⁻¹, and since 1 ml = 1 cm³ then we calculate c₀ = 6.7 x 10⁻⁶ M O₂ cm⁻³ [140]. Obviously this would be an overestimate, due to the fact that not all of the domain will be fully vascularised but at least we have obtained a reference value. The values of the non-dimensional parameters were given as

$$\begin{array}{l}{d}_n=0.0005,\text{a}\\ r=0.01,\text{b}\\ h=50,\text{c}\\ d_m=0.0005,\text{d}\\ k=1,\text{i}\\ s=0,\text{f}\\ d_c=0.5,\text{g}\\ n=0.5,\text{h}\\ w=0.57,\text{i}\\ f=\mathrm{0.025.}\text{j}\end{array}$$E34

7.4. Laser beam in a multiscale diffusion cancer-invasion model

7.4.1. Laser as a lorenz system

A laser system is the result of an interaction between the electromagnetic field and the substance, under certain circumstances. The Lorenz form of laser equations may be obtained using a semi-classical reasoning where the environment is analyzed quantically using the formalism of density matrix, and the electromagnetic field is treated classically, by means of Maxwell’s equations [141,142]. Here we consider only two energy levels of the involved microscopic systems (atoms, molecules, ions).

The first treatment of a two levels system was made by Bloch who analyzed the interaction of electrons with an oscillatory magnetic field superposed over a static magnetic field, in the framework of a magnetic resonance phenomenon. Due to the similarity of treatments and form of the obtained equations for the laser system, one can say that it forms the Maxwell-Bloch system.

Note that the density matrix method is applied in the treatment of laser systems, no matter how many energy levels, or number of oscillating modes are considered, as well as, in the consequent quantum treatment, where the electromagnetic field is quantized [143,144].

We start by discussing the effect of the electromagnetic field on the atoms of the environment. In the simplest situation, the electric field will induce in each atom an electric dipole whose moment is proportional to the field and is oriented along its direction. Neglecting the vectorial character, we have

$$\mu =\alpha E$$E35

where α is a constant characteristic to the type of the atom considered. If the concentration of (identical) atoms in the considered environment is N_a, then the polarization vector of the environment, equals the vectorial sum of all the dipole moments from the unit volume, and will be given by

$$P=N_a\mu ={\epsilon }_0\chi E$$E36

where ε₀ is the empty space permittivity and χ represents the electric susceptibility of the environment.

The problem of the induced dipole moment must be solved quantically using Schrödinger’s equation (see paragraph 5.2)

$$i\hslash \frac{\partial \psi }{\partial t}=\widehat{H}\psi$$E37

where $$\widehat{H}$$ is the Hamiltonian operator and ψ the wave function of the atom.

Since a monochromatic field of frequency ω₀, not very intense, interacts with the atom inducing transitions between two of its energetic levels, E₁ and E₂ i.e. E₂ - E₁ = ħω₀, it is usual to neglect the other levels and to approximate the atom as a system with two energy levels. If the wave functions of the atom in the two states are ψ₁ and ψ₂, respectively, then we have

$$\psi =C_1{\psi }_1+C_2{\psi }_2$$E38

where C₁, C₂ are the time dependent complex amplitude probabilities for the atom to find itself on the energy levels E₁ and E₂, respectively. In other words $${\rho }_{11}=C_1^{\ast }{C}_1\equiv {\left|C_1\right|}^2$$ represents the probability of the atom to find itself in the state ψ₁, and $${\rho }_{22}=C_2^{\ast }{C}_2\equiv {\left|C_2\right|}^2$$ the probability of the atom to find itself in the state ψ₂. The combinations $${\rho }_{12}=C_1{C}_2^{\ast }$$ and $${\rho }_{21}=C_2{C}_1^{\ast }$$ are transition probabilities between the two states. The four numbers ρ_ij (i, j = 1, 2) forms the so called density matrix for the 2-levels system. The asterisk attached to a parameter means the complex conjugate of the respective parameter. Obviously, ρ₂₁ = ρ^*₁₂.

The Hamiltonian operator of the system will consist of a sum between the Hamiltonian of the nonperturbed atom $${\widehat{H}}_0$$ and a term which describes the interaction of the atom with the field, $${\widehat{H}}^{\prime }$$,

$$\widehat{H}={\widehat{H}}_0+{\widehat{H}}^{\prime }$$E39

The functions ψ₁ and ψ₂ are eigenfunctions of the nonperturbed Hamiltonian, i.e.

$${\widehat{H}}_0{\psi }_1=E_1{\psi }_1,{\widehat{H}}_0{\psi }_2=E_2{\psi }_2,$$E40

Replacing (38), (39) and (40) into the Schrödinger equation, and after some standard calculus, we get the equations

$$\begin{array}{l}i\hslash {\dot{C}}_1=E_1{C}_1+{H^{\prime }}_{12}{C}_2,\text{a}\\ i\hslash {\dot{C}}_2=E_2{C}_2+{H^{\prime }}_{21}{C}_1,\text{b}\end{array}$$E41

where we introduced the notations

$$\begin{array}{l}{H^{\prime }}_{12}={\displaystyle \int {\psi }_1^{\ast }\widehat{H^{\prime }}{\psi }_2}dV,\text{a}\\ {H^{\prime }}_{21}={\displaystyle \int {\psi }_2^{\ast }\widehat{H^{\prime }}{\psi }_1}dV,\text{b}\end{array}$$E42

It has been taken into account the orthonormal property of the wave functions ψ₁ and ψ₂

$${\displaystyle \int {\psi }_i^{\ast }{\psi }_{j}dV={\delta }_{ij}}\left(i,j=1,2\right)$$E43

where δ_ij is Kroeneker’s symbol, and the fact that the interaction matrix $H_{ij}^{\text{'}}$ has no diagonal elements.

It is common to introduce the following simplification: if one chooses the zero energy value at the center of the interval between the two energies, then they become

$$\begin{array}{l}{E}_2=\left(1/2\right)\hslash {\omega }_0,\text{a}\\ E_1=-\left(1/2\right)\hslash {\omega }_0,\text{b}\end{array}$$E44

and Eqs. (41) rewrite

$$\begin{array}{l}i\hslash {\dot{C}}_1=-\frac{\hslash {\omega }_0{C}_1}{2}+{H^{\prime }}_{12}{C}_2\text{,}\text{a}\\ i\hslash {\dot{C}}_2=\frac{\hslash {\omega }_0{C}_2}{2}+{H^{\prime }}_{21}{C}_1,\text{b}\end{array}$$E45

In quantum mechanics, the electric dipole momentum is calculated as the expectation value of the classical electric dipole momentum μ = ex, where e is the electron charge and x is its displacement along the direction of the electric field.

For an atom in the state ψ this is given by

$$\mu ={\displaystyle \int {\psi }^{\ast }ex\psi dV}=\left(C_1^{\ast }{C}_2+C_1{C}_2^{\ast }\right){\displaystyle \int {\psi }_1^{\ast }ex{\psi }_{2}dV}$$E46

The integral $${\mu }_{12}={\displaystyle \int {\psi }_1^{\ast }ex{\psi }_{2}dV}$$ represents the electric dipole momentum of the interaction. In the considered approximation, the interaction Hamiltonian is identical to the classical expression of the interaction energy between an electric field and the induced electric dipole: U = -μE, but with μ $\to$ μ₁₂ (the dipole momentum of the interaction), i.e.

$${H^{\prime }}_{12}=-{\mu }_{12}E$$E47

By choosing a convenient phase relation between the wave functions, we can make μ₁₂ real so $H_{12}^{\text{'}}$, $H_{21}^{\text{'}}$ to be also real.

Eq. (46) suggests considering the expression

$$X=C_1{C}_2^{\ast }+C_1^{\ast }{C}_2\equiv {\rho }_{12}+{\rho }_{21}$$E48

We remark that the polarization (36) is expressed as a function of X through the equation

$$P=N_a{\mu }_{12}X=N_a{\mu }_{12}{\rho }_{12}+c.c.$$E49

where by c.c. we denote the complex conjugate of the previous expression.

We further consider the combinations

$$Y=i\left(C_1^{\ast }{C}_2-C_1{C}_2^{\ast }\right)\equiv -i\left({\rho }_{12}-{\rho }_{21}\right)$$E50

$$Z={\left|C_2\right|}^2-{\left|C_1\right|}^2\equiv {\rho }_{22}-{\rho }_{11}$$E51

For Z we also have a simple interpretation. If we multiply (51) by N_a we get the expressions N_aρ₂₂ and N_aρ₁₁, which represents the populations from the unit volume of the two levels, in other words, N_a Z ≡ N represents the difference of population between the levels (inversion of population) from the unit volume.

All the three expressions X, Y, Z are functions depending only on time. Their time derivatives are easily calculated using Eqs. (45) and their complex conjugates. The following relations result:

$$\dot{X}=-{\omega }_{0}Y-\frac{i}{\hslash }\left({H^{\prime }}_{12}-{H^{\prime }}_{21}\right)Z=-{\omega }_{0}Y$$E52

$$\dot{Y}={\omega }_{0}X+\frac{2{\mu }_{12}}{\hslash }EZ$$E53

$$\dot{Z}=-\frac{2{\mu }_{12}}{\hslash }EY$$E54

where in Eq. (52) the second form was obtained taking into account $H_{12}^{\text{'}}$ = $H_{21}^{\text{'}}$.

By multiplication of Eq. (52) with N_aμ₁₂, it transforms into an equation for $$\dot{P}$$, namely

$$\dot{P}=N_a{\mu }_{12}\dot{X}=N_a{\mu }_{12}\left({\dot{\rho }}_{12}+{\dot{\rho }}_{21}\right)=N_a{\mu }_{12}{\dot{\rho }}_{12}+c.c.$$E55

The equation for $${\dot{\rho }}_{12}$$ is obtained from Eqs. (52) and (53), taking into account Eqs. (48) and (50). It results

$${\dot{\rho }}_{12}=i{\omega }_0{\rho }_{12}+i\frac{{\mu }_{12}}{\hslash }EZ$$E56

Another equation for polarization is obtained by deriving Eq. (52) once again and using (53). It results

$$\ddot{P}+{\omega }_0^{2}P=-\frac{2{\omega }_0}{\hslash }{\mu }_{12}^{2}EN$$E57

It is interesting to note that, in the absence of the electromagnetic field, i.e. for E = 0, Eq. (57) describes a harmonic oscillator. This is unacceptable, since polarization is induced by the field, so it must attenuate after the field cancels. Physically, it occurs both because of the internal dynamics of the atomic (molecular) systems and of the dephasing between the oscillations of different dipoles by means of their self-interaction or their interaction with the crystal lattice (for a solid environment). This phenomenon is taken into consideration through phenomenological reasonings by introducing an amortization term of the form $$\dot{P}/T_2$$ in the equation. Eq. (57) transforms into the equation of a forced dumped oscillator:

$$\ddot{P}+\frac{\dot{P}}{T_2}+{\omega }_0^{2}P=-\frac{2{\omega }_0}{\hslash }{\mu }_{12}^{2}EN$$E58

Usually, the time T₂ is named transversal relaxation time. It is characteristic to the non-diagonal elements of the density matrix, so Eq. (56) must be rewritten

$${\dot{\rho }}_{12}+({\gamma }_{12}-i{\omega }_0){\rho }_{12}=i\frac{{\mu }_{12}}{\hslash }EZ$$E59

where γ₁₂ =1/T₂.

By multiplication of Eq. (54) with N_a the left side becomes $$\dot{N}$$. Replacing Y from (52) in (54), we get

$$\frac{\partial N}{\partial t}=\frac{2}{\hslash {\omega }_0}E\dot{P}$$E60

Eq. (60) shows that, at the disappearance of the electromagnetic field, the inversion of population must remain constant. However, an electromagnetic field resonant with the considered transition (ω ≈ ω₀) is composed of quanta which can be absorbed by atoms, so may have the effect of a transfer of population between the two levels. It is obvious that, at the canceling of the field, the inversion of population must evolve towards an equilibrium value N^e which is obtained by a process of pumping and by spontaneous relaxation processes. They imply the presence of other energetic levels besides those already considered. We proceed again by phenomenological reasonings. We suppose that this evolution is again exponentially, thus we add a term of the form $$(N-N^e)/T_1\equiv {\gamma }_{11}(N-N^e)$$, where T₁ is named longitudinal relaxation time (it is characteristic to the diagonal elements of the density matrix). The equation for the inversion of population becomes

$$\frac{\partial N}{\partial t}+{\gamma }_{11}(N-N^e)=\frac{2}{\hslash {\omega }_0}E\dot{P}$$E61

Eq. (61), together with (55) coupled with (59), or with (58) represents the substance equations. They must be associated with the electromagnetic field equation which we transcribe here

$${\nabla }^{2}E-\frac{{\eta }^2}{c^2}\frac{{\partial }^{2}E}{\partial t^2}={\mu }_0\frac{{\partial }^{2}P}{\partial t^2}$$E62

where c² = 1/ε₀μ₀, η =√ε/ε₀ is the refraction index of the environment (without the contribution of the transition between the two levels), P is the resonant part of the induced polarization and two non-conductive and non-magnetic laser environments were considered. In this case also, the energy losses produced by different mechanisms will be considered phenomenologically by means of an attenuation term introduced in the final form of the equation.

If we have a laser oscillator, the laser medium is placed between two mirrors which form an open cavity, and the oscillations may be triggered only by modes of oscillation characteristic to the cavity. They must satisfy the Helmholtz equation

$${\nabla }^{2}U+k^{2}U=0$$E63

where k = ηω_c/c, ω_c being the frequency of the considered mode (index c from cavity). For simplification, in what follows, we consider η=1 (gaseous environment).

We suppose the oscillation is produced on a single mode described by a spatial dependence of the form W(x, y, z). This dependence will characterize both the field and polarization, so we can write

$$E=\tilde{E}(t)W(x,y,z)\exp (i{\omega }_{c}t)+c.c.$$E64

and

$$P=\tilde{P}(t)W(x,y,z)\exp (i{\omega }_{c}t)+c.c.$$E65

respectively, where $$\tilde{E}(t)$$ and $$\tilde{P}(t)$$ are slowly time-varying (complex) amplitudes.

Since the vectors E and P have the same directions, we can neglect the vectorial aspect. Introducing (64) and (66) in (62), applying the slowly varying amplitude approximation, i.e. taking $$d^2\tilde{E}/d{t}^2\approx 0$$, $$d^2\tilde{P}/d{t}^2\approx -{\omega }_{c}P$$ and having in view that W(x, y, z) satisfies Eq. (63), we get

$$\frac{d\tilde{E}}{dt}=\frac{i{\omega }_c}{2{\epsilon }_0}\tilde{P}$$E66

It is necessary to include the loss by radiation which are due to, in the first place, mirrors imperfections. We do that by introducing a term $$\kappa \tilde{E}$$ in the left hand side, so the field equation becomes

$$\frac{d\tilde{E}}{dt}+\kappa \tilde{E}=\frac{i{\omega }_c}{2{\epsilon }_0}\tilde{P}$$E67

The equation for polarization is obtained comparing Eqs. (49) and (65). It results

$$\begin{array}{l}{\rho }_{12}=\frac{\tilde{P}W\exp (i{\omega }_{c}t)}{N_a{\mu }_{12}}\text{,}\text{a}\\ {\dot{\rho }}_{12}=\frac{1}{N_a{\mu }_{12}}(\frac{d\tilde{P}}{dt}+i{\omega }_c\tilde{P})W\exp (i{\omega }_{c}t)\text{b}\end{array}$$E68

which introduced in Eq. (59), and after multiplying both sides with W^*exp(-iω_ct) and integrating over the entire volume (mode) of the cavity, leads to

$$\frac{d\tilde{P}}{dt}+({\gamma }_{12}+i({\omega }_c-{\omega }_0))\tilde{P}=\frac{i{\mu }_{12}^2}{\hslash }\tilde{E}N$$E69

where, in the left hand side, we introduced the first term from (64) and where N represents the inversion of population from the volume occupied by the cavity mode, defined by the relation

$$N=\frac{N_a{\displaystyle \int ({\rho }_{22}-{\rho }_{11})W^{\ast }WdV}}{{\displaystyle \int W^{\ast }WdV}}$$E70

We then introduce the relations for the field and polarization (64) and (65) into the equation for the inversion (61). Using the approximation $$\dot{P}=i{\omega }_{c}P$$, neglecting the rapidly varying terms (which contain exp (±2iω_ct)), multiplying by WW^* and integrating over the volume of the cavity, we get

$$\frac{\partial N}{\partial t}+{\gamma }_{11}(N-N^e)=\frac{2i{A}_0}{\hslash }({\tilde{E}}^{\ast }\tilde{P}-\tilde{E}{\tilde{P}}^{\ast })$$E71

where N is given by Eq. (70). We consider ω_c/ω₀ = 1, and note

$$A_0=\frac{N_a{\displaystyle \int {(W^{\ast }W)}^{2}dV}}{{\displaystyle \int W^{\ast }WdV}}$$E72

Eqs. (67), (69) and (71) form the Bloch-Maxwell system and describe a unimodal laser oscillator. If we make the change of variables

$$\begin{array}{l}t=\frac{1}{{\gamma }_{12}}{t}^{\prime },\text{a}\\ \tilde{E}=\frac{i\hslash {\gamma }_{12}}{2\sqrt{A_0}{\mu }_{12}}A,\text{b}\\ \tilde{P}=\frac{\hslash {\gamma }_{12}{\epsilon }_0\kappa }{{\omega }_c\sqrt{A_0}{\mu }_{12}}R,\text{c}\\ N-N^e=\frac{2\hslash {\gamma }_{12}{\epsilon }_0\kappa }{{\omega }_c{\mu }_{12}^2}n\text{d}\end{array}$$E73

the Bloch-Maxwell equations gets the form of the Lorenz system. They become

$$\begin{array}{l}\frac{dA}{d{t}^{\prime }}=-\sigma A+\sigma R\text{a}\\ \frac{dR}{d{t}^{\prime }}=rA-R(1+i\Omega )-An\text{b}\\ \frac{dn}{d{t}^{\prime }}=-bn+\frac{1}{2}(A{R}^{\ast }+A^{\ast }R)\text{c}\end{array}$$E74

where the following notations were used

$$\begin{array}{l}\sigma =\frac{\kappa }{{\gamma }_{12}}\text{,}\text{a}\\ \Omega =\frac{{\omega }_c-{\omega }_0}{{\gamma }_{12}}\text{,}\text{b}\\ b=\frac{{\gamma }_{11}}{{\gamma }_{12}}\text{,}\text{c}\\ r=\frac{{\omega }_c{\mu }_{12}^2{N}^e}{2{\epsilon }_0\hslash \kappa {\gamma }_{12}}\text{d}\end{array}$$E75

In the form (74) the equations make up a complex Lorenz system. This was discussed in detail in the paper [145]. The complex Lorenz system transforms into the well known real Lorenz system if the resonance is exact (ω_c = ω₀ ⇒ Ω = 0) and the phases of the amplitudes $$\tilde{E}$$ and $$\tilde{P}$$ are chosen so the functions A and R to be real

$$\begin{array}{l}\frac{dA}{d{t}^{\prime }}=\sigma (R-A)\text{a}\\ \frac{dR}{d{t}^{\prime }}=A(r-n)-R\text{b}\\ \frac{dn}{d{t}^{\prime }}=AR-bn\text{c}\end{array}$$E76

The fact that a unimodal laser oscillator is described by the Lorenz system was remarked for the first time by Haken (Haken, 1975). Therefore, it was demonstrated that the immense variety of dynamical behaviors, including the chaotic ones, presented by a Lorenz system, must be expected to occur in a laser. Among the first who reported Lorenz type chaotic behaviors in a laser, were Weiss and Brock [147].

Equations of the same form are obtained also when one considers a travelling wave laser, such as laser amplificators where the wave passes only one time the environment, or a circular unidirectional laser [148,149].

7.5. A chaotic multi-scale cancer-invasion model

From the non–dimensional space-time model (33), discretization was performed by neglecting all the spatial derivatives resulting in the following simple 4D temporal dynamical system

$$\begin{array}{l}\frac{dn}{dt}=0\text{a}\\ \frac{df}{dt}=-\eta mf\text{b}\\ \frac{dm}{dt}=\kappa n-\sigma m\text{c}\\ \frac{dc}{dt}=\nu f-\omega n-\phi cd\end{array}$$E77

When simulated, the temporal system (77) with the set of parameters (34) exhibits a virtually linear temporal behavior with almost no coupling between the four concentrations that have very different quantitative values (all phase plots between the four concentrations, not shown here, are virtually one-dimensional). To see if a modified version of the system (77) could lead to a chaotic description of tumor growth, and following the method in [150], four new parameters, a₁, a₂, a₃, and a₄ are introduced. The resulting model is

$$\begin{array}{l}\frac{dn}{dt}=0\text{a}\\ \frac{df}{dt}=a_1\eta (m-f)\text{b}\\ \frac{dm}{dt}=f\left(a_3-c\right)-m+a_2\kappa n\text{c}\\ \frac{dc}{dt}=\nu fm-a_4\phi c-\omega n\text{d}\end{array}$$E78

The introduction of the parameters (a₁, a₂, a₃, a₄) was motivated by the fact that tumor cell shape represents a visual manifestation of an underlying balance of forces and chemical reactions [151]. Specifically, the parameters represent the following quantities: a₁ = tumor cell volume (proliferation/non-proliferation fraction), a₂ = glucose level, a₃ = number of tumor cells, a₄ = diffusion from the surface (saturation level).

A tumor is composed of proliferating (P) and quiescent (or non-proliferating) (Q) cells. Tumor cells shift from class P to class Q as the tumor grows in size [152]. Model dependence on the ratio of proliferation to non-proliferation is introduced via the first parameter, a₁. The discretization of Eq. (33a) leads to cell density being modeled as a constant in Eq. (78a). Accordingly, cell density does not play a role in the dynamics. In (78) the cell density is re-introduced into the dynamics via the cell number, a₃. The importance of introducing a₃ also appears in connection with the cyclin-dependent kinase (Cdk) inhibitor p27, the level and activity of which increase in response to cell density. Levels and activity of Cdk inhibitor p27 also increase with differentiation following loss of adhesion to the ECM [153].

The ability to estimate the growth pattern of an individual tumor cell type on the basis of morphological measurements should have general applicability in cellular investigations, cell–growth kinetics, cell transformation and morphogenesis [154].

Cell spreading alone is conducive to proliferation and increases in DNA synthesis, indicating that cell morphology is a critical determinant of cell function, at least in the presence of optimal growth factors and extracellular matrix (ECM) binding [155]. The varying morphology of most cells can stimulate cell proliferation through integrin-mediated signaling indicating that cell shape may govern how individual cells will respond to chemical signals [156].

Parameters (a₁, a₂, a₃, a₄), introduced in connection with cancer cells morphology and dynamics could also influence the very important factor chromatin associated with aggressive tumor phenotype and shorter patient survival time.

For computations, the parameters were set to a₁=0.06, a₂=0.05, a₃=26.5 and a₄=40. Small variation of these chosen values would not affect the qualitative behavior of the new temporal model (78). Simulations of (78), using the same initial conditions and the same non-dimensional parameters as before, show chaotic behavior in the form of Lorenz-like strange attractor in the 3D (f−m−c) subspace of the full 4D (n−f−m−c) phase-space (Figs. 13-Figs. 16).

Figure 13.

A 3D Lorenz-like chaotic attractor from the modified tumor growth model (78 b-d). The attractor effectively couples the MM–concentration f, the MDE–concentration m, and the oxygen concentration c in a mask–like fashion.

Figure 14.

The m − f phase plot of the 3D attractor.

Figure 15.

The m − c phase plot of the 3D attractor.

Figure 16.

The f − c phase plot of the 3D attractor.

The space-time system of rate PDEs corresponding to the system in (33) provides the following multi-scale cancer invasion model

$$\begin{array}{l}\frac{\partial n}{\partial t}=d_n{\nabla }^{2}n-\rho \nabla \cdot \left(n\nabla f\right)\text{a}\\ \frac{\partial f}{\partial t}=a_1\eta (m-f)\text{b}\\ \frac{\partial m}{\partial t}=d_m{\nabla }^{2}m+a_{2}kn+f\left(a_3-c\right)-m\text{c}\\ \frac{\partial c}{\partial t}=d_c{\nabla }^{2}c+\nu fm-\omega n-a_4\phi c\text{d}\end{array}$$E79

The new tumor–growth model (79) retains all the qualities of the original model (33) plus includes the temporal chaotic ‘butterfly’–attractor. This chaotic behavior may be a more realistic view on the tumor growth, including stochastic–like long–term unpredictability and uncontrollability, as well as sensitive dependence of a tumor growth on its initial conditions.

Now, if we compare (78) with (76) and make a one-to-one correspondence between these two systems of equations, we see that A which is the electric field amplitude corresponds to f, the MM concentration, R which is the polarization amplitude corresponds to m, the MDE concentration and n the inversion of population corresponds to c, the oxygen concentration. Since both systems, the laser and the tumor invasion can be written in the form of a Lorenz system, we can suppose that the metastatic cancer cells moving through the systemic circulation form a coherent wave, i.e. a particular type of chemically pumped (since it may obtain its energy from chemical reactions) laser with oxygen. In the following section we show moreover, that this coherent wave can be identified with a travelling wave laser with oxygen.

7.6. Travelling waves in the multiscale diffusion cancer-invasion model

Let us write the system (28) again. We assume the model refers to the the averaged behavior of the tumor cells in the direction of invasion only and ignores variations in a plane normal to the direction of invasion.

Invasive cells. Since in their experiments Aznavoorian et al. [157] reported minimal chemokinetic movement, a key feature of the following model is the absence of the term for random cell motility. Also, we introduce a term of increased proliferation of malignant cells relative to normal cells, F(n), which will be initially modeled as a logistic type growth of the form k₁ n (k₂ − n) which has been shown [158], in order for us to describe adequately the growth of human tumors grown [159].

Extracellular matrix. The motility of extracellular matrix elements, unlike the one of malignant cells and oxygen, is negligible due to the fact that these elements are much longer than cells. The dynamics of connective tissue can therefore be modeled as a simple passive degradation by the activity of the tissue proteases; this proteolysis can now be described by −G(f, m), since it depends on the amount of collagen f still present as well as on the protease m.

Proteases. Proteases generation is narrowly confined to the interface invading tumor and receding connective tissue interface. In some cases it is possible to localize the interstitial collagenase production to the stromal fibroblasts immediately adjacent to the site of tumor invasion, possibly leading to the fact that invasive cells release a stimulus for induction of interstitial collagenase by fibroblasts. Nabeshima et al. [160] managed to sequence a tumor cell derived collagenase stimulatory factor. Protease generation located only at the invading front can be explained however in other ways. In their work, Xie et al. [161] have revealed that the induction of 92-kd type IV collagenase activity in cultures of A431 human epidermoid carcinoma cells is density dependent. They showed that only dividing cells stained positive when treated with anti-MMP antibodies and as a consequence only noncontact-inhibited tumor cells produce protease. Many proteases are predominantly membrane bound (e.g. uroplasminogen activator), but even when the protease is secreted into the extracellular space, activation occurs only on the cell surface, so as a result the behavior closely resembles that for membrane bound proteases [162]. Therefore protease diffusion in the model is not included. We must then define the function H(n, f) to represent the dependence of this tightly regulated protease production on the local concentrations of the melanoma cells and collagen. In addition we assume that the protease decays linearly, with half-life K.

Oxygen. As in the original model, we presuppose it diffuses into the MM, decays naturally, is consumed by the tumor and for simplicity, oxygen production is proportional to the MM density. Therefore, we introduce the function I(n, f) and c decays linearly, with half-life Λ. The parameter c does not appear anywhere else in the system, so this equation will be easily separated.

Combining all of the above, we are now ready to write the model as

$$\begin{array}{l}\frac{\partial n}{\partial t}=\stackrel{invasive\text{ cell proliferation}}{\overbrace{F(n)}}-\stackrel{haptotacticcell\text{ movement}}{\overbrace{k_3\frac{\partial }{\partial x}\left(n\frac{\partial f}{\partial x}\right)}}\text{a}\\ \frac{\partial f}{\partial t}=-\stackrel{proteolysis}{\overbrace{G(f,m)}}\text{b}\\ \frac{\partial m}{\partial t}=\stackrel{proteaseproduction}{\overbrace{H(n,f)}}-\stackrel{naturaldecay}{\overbrace{Km}}c\\ \frac{\partial c}{\partial t}=\stackrel{diffusion}{\overbrace{d_c\frac{{\partial }^{2}c}{\partial x^2}}}+\stackrel{productionanduptake}{\overbrace{I(f,n)}}-\stackrel{decay}{\overbrace{\Lambda c}}\text{d}\end{array}$$E80

where F, G, H and I are functions of n, f and m. Compared to previous work on the modeling of cell movement, this model is unusual in that there is no cellular diffusion. This case has been considered previously [163] in the very different context of cellular aggregation, where they obtained conditions for blow-up in the absence of cell kinetics.

Before proceeding further, we must eliminate m from the equations as follows. The time scales associated with protease production and protease decay are much shorter than a typical timescale for the invading cells. Hence writing $$H(n,f)=K\overline{H}(n,f)$$, where we assume K>>1, and multiplying through Eq. (80c) by the small parameter K⁻¹, we deduce that to leading order $$m=\overline{H}(n,f)$$. Henceforth no reference to m is needed: this expression may be used to eliminate m from Eqs. (80a) and (80b). In the same way, writing $$I(f,n)=\Lambda \overline{I}(f,n)$$, assuming Λ >> 1 and multiplying through Eq. (80d) by the small parameter Λ⁻¹, we deduce that to leading order $$c=\overline{I}(f,n)$$. This type of quasi-steady state assumption is a common one in enzyme kinetics [164], and numerical simulations of the four equations (80a-d) are in good accordance with the simplified system of two equations; a strong point of the two equation case is that it is amenable to detailed mathematical analysis.

We examine the model using the simple functional forms

$$\begin{array}{l}F(n)=k_{1}n(k_2-n),\text{a}\\ G(f,m)=k_{4}mf,\text{b}\\ H(n,f)=k_{5}nf,\text{c}\\ I(f,n)=k_{6}f-k_{7}n\text{d}\end{array}$$E81

7.7. Nondimensionalization

After making the substitutions for F, G, H and I from (81) into equations (80a-d) and eliminating m using m = H(n, f) and c using c = I(f, n) we nondimensionalize the resulting equations using

$$\tilde{n}=\frac{n}{n^{\ast }},\text{ <mover accent="true">f <mo>˜</mo></mover>}=\frac{f}{f^{\ast }},\text{ <mover accent="true">t <mo>˜</mo></mover>}=\frac{t}{T},\text{ <mover accent="true">x <mo>˜</mo></mover>}=\frac{x}{L},L={\left[\frac{k_3}{k_2{k}_4{k}_5}\right]}^{1/2},T=\frac{1}{k_1{k}_2},n^{\ast }=k_2,f^{\ast }=\frac{k_1}{k_4{k}_5}$$

Dropping tildes for notational convenience then gives rise to the system

$$\begin{array}{l}\frac{\partial n}{\partial t}=n(1-n)-\frac{\partial }{\partial x}\left[n\frac{\partial f}{\partial x}\right]\text{a}\\ \frac{\partial f}{\partial t}=-n{f}^2\text{b}\end{array}$$E82

7.8. Spatially homogeneous system

Setting ∂/∂x = 0 in (82a) we note that the spatially homogeneous system has two steady states:

1. n = 0, f arbitrary − this is a continuum of (unstable) steady states parameterized by the (variable) amount of connective tissue in different tissues;

2. n = 1, f = 0 − this (stable) steady state corresponds to complete replacement of the normal tissue by invading malignant cells.

With ∂/∂x = 0, (82a) and (82b) can be solved explicitly giving

$$\begin{array}{l}n(t)={\left[1+\exp (-t+c_2)\right]}^{-1}\text{a}\\ f(t)={[}^t\text{b}\end{array}$$E83

where c₁ and c₂ are arbitrary constants. The behavior as t $\to$ ∞ shows that n $\to$ 1 and f $\to$ 0, hence justifying our classification of the steady state (n, f) = (1, 0) as stable.

7.9. Travelling wave analysis

The invasion process should normally correspond to the travelling wave solutions of the model (82a) and (82b) with the normal tissue steady state n = 0 ahead of the wave and the fully malignant state n = 1, f = 0 behind the wave. This fact is verified by numerical solutions of (82a) and (82b), which are not detailed here. These travelling wave solutions can be studied analytically using the travelling wave differential equations. We look for constant shape travelling wave front solutions of (82a) and (82b) by setting

$$\begin{array}{l}n(x,t)=N(z),\text{a}\\ f(x,t)=F(z),\text{b}\\ z=x-\xi t\text{c}\end{array}$$E84

where ξ is the positive wave speed which has to be determined. When solutions of the type (84) exist, they represent travelling waves moving in the positive x - direction. Substitution of (84) into (82a) and (82b) followed by simple algebraic manipulation gives

$$\begin{array}{l}\frac{dN}{dz}\left(-\xi +\frac{2N{F}^2}{\xi }\right)=N(1-N)-\frac{2N^3{F}^3}{{\xi }^2}\text{a}\\ \frac{dF}{dz}=\frac{N{F}^2}{\xi }\text{b}\end{array}$$E85

The analysis of (85a,b) involves the study of the (N, F) phase plane. Since we are looking for travelling waves connecting (1, 0) and $$(0,\widehat{F})$$ in the (N, F) phase plane we look for solutions of (85a,b) with boundary conditions

$$\begin{array}{l}N(-\infty )=1,\text{a}\\ F(-\infty )=0,\text{b}\\ N(\infty )=0,\text{c}\\ F(\infty )=\widehat{F},\text{d}\end{array}$$E86

which requires (1, 0) to have an unstable manifold while $$(0,\widehat{F})$$ must have a stable manifold. In order to study this we look at the stability of the system (85a,b).

7.10. Stability analysis

The steady states (N₀, F₀) of (85a) and (85b) are $$(0,\widehat{F})$$ and (1, 0), where $$\widehat{F}$$ represents a steady states continuum. We study their stability by analyzing the eigenvalues of the stability matrix linearized about the steady states. The eigenvalues about $$(0,\widehat{F})$$ are −1/ξ and 0. The corresponding eigenvectors are $$(1,-\widehat{F})$$ and (0, 1). The negative eigenvalue shows that there is a stable manifold along $$(1,-\widehat{F})$$. The zero eigenvalue represents translations along the steady states continuum.

The eigenvalues about (1, 0) are 1/ξ and 0. The eigenvector corresponding to 1/ξ is (1, 0) and represents movement along the N axis. The eigenvector corresponding to the zero eigenvalue is (0, 1) which is in the direction normal to the N axis. The trajectory leaving this steady state leaves along the eigenvector corresponding to the zero eigenvalue of the linearized system, as shown by the numerical solutions of (85a,b). This zero eigenvalue is a result of (85b). In order to get a more detailed and clear image of the behavior close to (1, 0) we must look at the nonlinear terms in (85b). One way to do this is to use the techniques of the centre manifold theory which shows that as z $\to$ - ∞, N(z) approaches 1 exponentially while F(z) tends zero as z⁻¹.

The existence of an unstable manifold about $$(0,\widehat{F})$$ as z $\to$ ∞ and a stable centre manifold about (1, 0) as z $\to$ −∞ is consistent with the existence of a travelling wave orbit connecting the two steady states.

8. Conclusions

1. Cancer cannot be reduced to simple mathematical principles. Its irregular mode of carcinogenesis, erratic tumor growth, variable response to tumoricidal agents, and less-known metastatic patterns constitute highly variable clinical behavior. Characterizing this process requires an accurate understanding of tumor cells and host tissues interactions and ultimately determines prognosis. Applying time-tested and evolving mathematical methods to oncology may provide new methods, with inherent advantages, for the description of tumor behavior, selection of therapeutic modes, prediction of metastatic patterns, and the defining of an inclusive basis for prognostication. Mathematicians describe equations that define tumor growth and behavior, whereas surgeons actively deal with biological processes. Mathematics in oncology applies these principles to clinical settings.

2. The main conclusions of this work are as follows:

   1. mathematics of cancer proves to be chaotic and highly non-linear, justifying the use of space(-time) non-differentiability as a starting base model;

   2. a chaotic multi-scale cancer-invasion model is manufactured, which embeds a Lorenz attractor in its solutions;

   3. since laser can be expressed as a Lorenz system, we may assume some correspondences between the laser and the above mentioned chaotic multi-scale cancer-invasion model;

   4. the basic model for solid tumor growth admits a travelling wave solution;

   5. we suggest that metastatic tumor cells which move through the systemic circulation are similar to a coherent wave, i.e. a travelling wave chemically pumped type laser with oxygen;

   6. we assume the extracellular matrix and in particular, the tumor microenvironment are non-differential media endowed with holographic properties (capacity to memorize, interference abilities and source of forces);

   7. the two well-known phenomena: tumor self-seeding by CTC and hypoxia, in our opinion, both support the idea of complete holography (a hologram which becomes the very object in the particular case of living organisms).

1. Experimentally testable, mathematics applied in oncology may provide a framework to determine clinical outcome on a patient-specific basis and increase the growing awareness that mathematical models help simplify seemingly complex and random tumor behavior.