AN IMAGING-DRIVEN, MECHANICAL DEFORMATION-COUPLED REACTION-DIFFUSION MODEL FOR DESCRIBING TUMOR DEVELOPMENT

Compressive stresses play important roles on tumor cells proliferation and invasion. In this work we propose to inform a tumor growth model with the recovered deformation coming from in vivo imaging data. We investigated the use of an optical flow technique, the well known Lucas-Kanade method, to capture tumor deformation in a synthetic experimental breast cancer setting. We compare displacements and stresses obtained with this method with those derived from a previously developed reaction–diffusion model with mechanical deformation. We show that the considered optical flow technique may capture deformations appearing in breast cancers, being a useful alternative to integrate in vivo data to mathematical tumor models. Palavras-chave: Mathematical model. Tumor growth. Deformation. MRI. Optical flow. Resumo: Tensões de compressão desempenham papéis importantes na proliferação e invasão de células tumorais. Neste trabalho propomos informar um modelo de crescimento tumoral com a deformação proveniente de dados de imagem in vivo. Investigamos o uso de uma técnica de fluxo ótico para capturar deformações tumorais em um cenário sintético de câncer de mama. Comparamos os deslocamentos e tensões obtidas com esta técnica com aquelas derivadas de um modelo de reação-difusão com deformação mecânica. Mostramos que a técnica de fluxo ótico considerada pode capturar deformações que aparecem em cânceres de mama, sendo uma alternativa útil para integrar dados in vivo a modelos que visam representar o crescimento de tumores. Palavras-chave: Modelo matemático. Crescimento tumoral. Deformação, MRI, Fluxo ótico. 1M.Sc. Computational Modeling, Laboratório Nacional de Computação Científica, annacmr@lncc.br 2D.Sc. Computational Modeling, Universidade Federal de Juiz de Fora, rafael.bonfim@ice.ufjf.br 3D.Sc. Computational Modeling, Oden Institute for Computational Eng. and Sciences, lima@ices.utexas.edu 4M.Sc. Computer Science, Laboratório Nacional de Computação Científica, naozuka@lncc.br 5D.Sc. Nuclear Engineering, Laboratório Nacional de Computação Científica, rcca@lncc.br Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 217-01, 217-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21142 1


INTRODUCTION
As a solid tumor evolves, compressive stresses accumulate within the tumor due to growth. These stresses play important roles on tumor cells phenotype differentiation and tumor microenvironment conditions (CHENG et al., 2009;STYLIANOPOULOS et al., 2012). Specifically, higher compressive stresses are known to inhibit cell proliferation and stimulate invasion and apoptosis; to decrease perfusion due to compression of blood vessels, and to increase hypoxia (STYLIANOPOULOS et al., 2012). Many mathematical models have been developed in the literature to represent tumor growth under deformation (WEIS et al., 2013;LIMA et al., 2016;VOUTOURI;STYLIANOPOULOS, 2018). From a continuum mechanics point of view, they are usually built by performing a kinematic decomposition, applying the momentum balance equation, and adding constitutive relations. This framework involves a series of assumptions that ultimately impact the prediction of the tumor deformation.
A different framework can be pursued by using in vivo data to recover the tumor deformation. It is now recognized that "advanced quantitative imaging measurements from [magnetic resonance imaging] MRI or positron emission tomography can provide quantification of various tumor characteristics" (JARRETT et al., 2018). Moreover, non-invasive imaging measurements allow identifying patient-specific characteristics which would provide a more feasible picture of individualized tumor growth and response. Thus, here we investigate the use of an optical flow technique to track tumor deformation. The main idea is to evaluate at each image pixel the stress tensor which will be used to feed, in the future, a desired tumor growth model through a stress-dependent proliferation or migration feature. Basically, given two consecutive images where the tumor is observed, we first define the image brightness associated with each image pixel and time. Assuming brightness conservation and piecewise smooth motion, there are many methods to recover the optical flow that allows determining the deformation tensor. Here we investigate the Lucas-Kanade (LUCAS;KANADE, 1981) technique. We then compare the proposed approach to that one based on a previously developed tumor model built on the continuum mechanics framework (LIMA et al., 2016), for a synthetic experimental breast cancer setting. We show that optical flow techniques may capture deformations appearing in breast cancers, being a useful alternative to integrate in vivo data to mathematical tumor models.
The manuscript is outlined as follows. In Section 2 we present a mechanical deformationcoupled reaction-diffusion model for describing tumor development. Section 3 introduces the Lucas-Kanade method. In Section 4 we define a synthetic breast cancer scenario to be used to compare the continuum mechanics and optical flow approaches. The last section summarizes our conclusions. Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 217-01, 217-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21142

MATHEMATICAL MODEL
Here we use the reaction-diffusion model with elastic deformation investigated in (LIMA et al., 2016), in which a small displacement regime is assumed. Defining the displacement field by u, the deformation (strain) tensor is given by: where F is the deformation gradient tensor, and I the identity tensor. The deformation tensor can be decomposed into the elastic and growth parts as: in which ε s represents the mechanical counterpart of strain due to applied stress and ε g stands for the stress-free strain due to growth. We also define the linear elastic inhomogeneous material tensor as: where φ stands for the tumor cell density and g(φ) is an interpolation function that allows C(φ) to vary smoothly across the interface from the healthy tissue (C h ) to the tumor tissue WISE et al., 2008). As a remark, (WEIS et al., 2013) state that the breast tissue is often modeled using a viscoelastic assumption. For this preliminary study on the tumor mechanical effects, we assume a linear elastic material model. We also do not distinguish between the mechanical properties of the tumor microenvironment. The model development continues by defining the strain energy density as: in which σ(φ) = λφI is a linear isotropic compositional stress tensor associated with the tumor growth, with λ depending on the tumor proliferation. Considering a quasi-static approximation for the linear momentum balance and assuming no body forces, Eq. (4) yields: Rewriting Eq. (5) in terms of the displacement vector u, we obtain: Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 217-01, 217-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21142 Finally, for an elastic, isotropic and homogeneous tumor microenvironment, and defining the Poisson's ratio ν, the Young's modulus E, and the shear modulus G = E/(2(1 + ν)), the linear momentum equilibrium Eq. (6) becomes: The tumor cell dynamics is modeled by a reaction-diffusion equation. By assuming that tumor cells diffuse randomly in the microenvironment and grow bounded by the available resources, the cell balance equation is given by: where D is the tumor cell diffusion coefficient and k is the tumor cell maximum rate of proliferation. The coupling between models (7) and (8) is performed by a negative feedback control on the tumor proliferative growth rate. Specifically, we assume that the maximum rate of proliferation exponentially decays at a rate depending on the von Mises stress (WEIS et al., 2013), which yields: where γ is a scaling constant, σ vm the von Mises stress, andk is the maximum proliferation rate in a free from stresses scenario.
To summarize, the tumor mathematical model with mechanical coupling is defined on Ω × (0, τ max ] and is given by: with all parameters defined in Table 1 for the proposed two-dimensional synthetic breast cancer scenario. Implementation details as well as information on Ω, τ max , initial and boundary conditions are given in Section 4. We remark that identifying in vivo parameters is important to allow new therapies to be developed, as well as to provide insights into the mechanisms of tumor growth (KATIRA; BONNECAZE; ZAMAN, 2013). In particular, our aim in this work is to automate the acquisition of tumor deformation through the use of the optical flow technique. Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 217-01, 217-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21142

OPTICAL FLOW
The automatic quantification of physiological parameters (such as strain rate, stress state, mechanical properties of tissues, among others) from image exams can be done through the optical flow technique. By definition, optical flow is the motion of the brightness patterns in a sequence of images (HORN; SCHUNCK, 1981). This technique was introduced by (HORN; SCHUNCK, 1981) and (LUCAS; KANADE, 1981), and many variations have been developed in the literature since then. Here, we investigate the use of the Lucas-Kanade approach to represent the stress state in a breast cancer scenario. To introduce the method, assume that I(x, y, t) is a function describing the time sequence of two-dimensional images, where I represents the intensity of the pixel located at position x = (x, y) of the image at time t. The Lucas-Kanade method determines the velocity v of each pixel of the image by minimizing the functional: where ∇I = (I x , I y ) is the spatial intensity gradient and W(x) denotes a 2D Gaussian mask that emphasizes more influence to constraints at the center of the neighborhood than those at the boundary.
The Lucas-Kanade method assumes that the velocity v = (u, v) of each pixel of the image is constant within a spatial neighborhood B of the pixel under consideration and the material derivative DI/Dt vanishes over the entire time sequence.
For n pixels in the positions x i ∈ B at a single time t, the minimization of (11) results in the linear system: in which The solution of (12) is as shown by (BARRON; FLEET; BEAUCHEMIN, 1994). We adopted n = 5 neighborhoods in the algorithm implemented in this work. The partial derivatives I x , I y and I t are approximated according to (HORN;SCHUNCK, 1981).

NUMERICAL EXPERIMENTS
In the following, we use a synthetic experimental breast cancer setting to investigate tumor growth and invasion through the two-dimensional coupled model given in Section 2. The approximate solution to system (10) is solved by using the finite element method and the implicit Euler method to approximate the time derivatives. The system is linearized by using Picard's method and is solved uncoupling the equations through a Gauss-Seidel strategy. The resulting discrete forms are implemented in C ++ using the open source libMesh library (KIRK et al., 2006). We also define Ω and τ max to complete the definition of our spatial and temporal domains, respectively. An example of MRI can be seen in Fig. 1(a) for a specific patient (HUANG et al., 2014a;HUANG et al., 2014b;CLARK et al., 2013). Figure 1(b) shows a synthetic breast domain discretization partitioned into 6480 elements. The temporal domain (0, τ max ] is split into uniform time steps ∆t = 0.1 d, with τ max = 10 d. Our discretization settings are given in Table 2. The initial condition for tumor cells is assumed to be of the following circular shape:   (HUANG et al., 2014a;HUANG et al., 2014b;CLARK et al., 2013) and the breast domain discretization used in this work.
The dynamics is completely described by defining the following homogeneous Neumann and Dirichlet conditions: where n is a unit exterior normal vector on ∂Ω.
The tumor evolution, the magnitude of the displacement and the von Mises stress are shown in Fig. 2 for two different time points, t = 1.0 d (first row) and t = 1.5 d (second row). Figure 2(a) depicts the tumor concentration in the initial phase of the tumor evolution, whereas Fig. 2(d) shows how the tumor profile evolved 12 hours later. As expected, tumor growth is isotropic. Figures 2(b) and 2(e) show the magnitude of the displacement as a measure of how the tumor growth experienced in both time points caused deformations to the tumor's surrounding tissue. We also remark the von Mises stress in Figs. 2(c) and 2(f) to reflect the total stress experienced in the breast tissue.
In order to investigate how the Lucas-Kanade method behaves we compare the displacements and stresses previously obtained with the ones calculated by the technique described in Section 3. In the proposed approach, tumor motion is tracked using Figs. 2(a) and 2(d). With both images and the definition of the brightness associated with each pixel and time we can recover the displacement suffered from one image to another to then calculate information such as the von Mises stress. Figures 3(b) and 3(d) show the magnitude of the Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 217-01, 217-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21142 displacement and the von Mises stress calculated with the displacement obtained by the Lucas-Kanade method. Since Lucas-Kanade captures changes in the pixels of the considered images, it may not detect regions whose tumor concentration is too low, therefore we define an interest region to evaluate our model solutions. We stated that it is only appropriate to look at regions where tumor concentration is above a certain threshold, as e.g., ε = 1×10 −2 . With this value we  Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 217-01, 217-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21142 calculated the difference between the displacements and the von Mises stresses in the selected region. These results are shown in Figs. 3(a) and 3(c). Comparing both model and optical flow outcomes we can see that the results are similar in magnitude, although there are some discrepancies. Of note, the slight lack of symmetry in the tumor growth was not completely captured by the present optical flow technique. To better quantify those differences, we consider different ε and calculate the absolute error between the values provided by each method. The results can be seen in Fig. 4 for ε = 1 × 10 −2 , ε = 1 × 10 −3 and ε = 1 × 10 −4 . The analysis points out that using ε = 1 × 10 −2 the Lucas-Kanade method better captures the deformations suffered in this breast cancer setting. The difference between the results seems to occur due to small tumor concentration fluctuations in the selected region which yield small differences in the brightness of the image that are not precisely captured by the optical flow method. Notice that the optical flow method takes into account the difference between pixels of two images (discrete approach), which leaves out of the equation those slight variations in concentration.
In contrast, the displacements in the model are still calculated even in the regions where tumor concentrations are very low (see Fig. 2(b) and Fig. 2(e) for reference). The definition of ε controls this feature to a certain extent, but it is a limited alternative. This is a point worth investigating in a future work. Other current optical flow methods like Farnebäck, TV-L1, and EpicFlow (ALLAERT et al., 2019) will also be examined.

CONCLUSIONS
In this work we show that the Lucas-Kanade method may capture deformations appearing in breast cancers, being a useful alternative to integrate in vivo data to mathematical tumor models. Such approach is particularly useful since it allows integrating patient-specific image information to cancer models, ultimately providing a personalization of the model predictions. Besides investigating other optical flow techniques, next step in the present research will pursue to inform our models with the recovered deformation coming from in vivo dynamic contrastenhanced MRI data at different time steps. In this way we will be able to include in our models patient-specific characteristics which would provide a more feasible picture of individualized tumor growth and response to treatment.   The results are displayed for ε = 1 × 10 −2 , ε = 1 × 10 −3 and ε = 1 × 10 −4 .