A dynamic model of CT scans for quantifying doubling time of ground glass opacities using histogram analysis

  • Received: 09 October 2017 Revised: 20 January 2018 Published: 01 October 2018
  • MSC : Primary: 92C50, 92C55; Secondary: 35K58

  • We quantify a recent five-category CT histogram based classification of ground glass opacities using a dynamic mathematical model for the spatial-temporal evolution of malignant nodules. Our mathematical model takes the form of a spatially structured partial differential equation with a logistic crowding term. We present the results of extensive simulations and validate our model using patient data obtained from clinical CT images from patients with benign and malignant lesions.

    Citation: József Z. Farkas, Gary T. Smith, Glenn F. Webb. A dynamic model of CT scans for quantifying doubling time of ground glass opacities using histogram analysis[J]. Mathematical Biosciences and Engineering, 2018, 15(5): 1203-1224. doi: 10.3934/mbe.2018055

    Related Papers:

    [1] Yuan-Nan Young, Doron Levy . Registration-Based Morphing of Active Contours for Segmentation of CT Scans. Mathematical Biosciences and Engineering, 2005, 2(1): 79-96. doi: 10.3934/mbe.2005.2.79
    [2] Anam Mehmood, Ishtiaq Rasool Khan, Hassan Dawood, Hussain Dawood . A non-uniform quantization scheme for visualization of CT images. Mathematical Biosciences and Engineering, 2021, 18(4): 4311-4326. doi: 10.3934/mbe.2021216
    [3] Javad Hassannataj Joloudari, Faezeh Azizi, Issa Nodehi, Mohammad Ali Nematollahi, Fateme Kamrannejhad, Edris Hassannatajjeloudari, Roohallah Alizadehsani, Sheikh Mohammed Shariful Islam . Developing a Deep Neural Network model for COVID-19 diagnosis based on CT scan images. Mathematical Biosciences and Engineering, 2023, 20(9): 16236-16258. doi: 10.3934/mbe.2023725
    [4] Xu Bie, Yuanyuan Tang, Ming Zhao, Yingxi Liu, Shen Yu, Dong Sun, Jing Liu, Ying Wang, Jianing Zhang, Xiuzhen Sun . Pilot study of pressure-flow properties in a numerical model of the middle ear. Mathematical Biosciences and Engineering, 2020, 17(3): 2418-2431. doi: 10.3934/mbe.2020131
    [5] Vasileios E. Papageorgiou, Georgios Petmezas, Pantelis Dogoulis, Maxime Cordy, Nicos Maglaveras . Uncertainty CNNs: A path to enhanced medical image classification performance. Mathematical Biosciences and Engineering, 2025, 22(3): 528-553. doi: 10.3934/mbe.2025020
    [6] Huiying Li, Yizhuang Song . Sparse-view X-ray CT based on a box-constrained nonlinear weighted anisotropic TV regularization. Mathematical Biosciences and Engineering, 2024, 21(4): 5047-5067. doi: 10.3934/mbe.2024223
    [7] Dayu Xiao, Jianhua Li, Ruotong Zhao, Shouliang Qi, Yan Kang . Iterative CT reconstruction based on ADMM using shearlet sparse regularization. Mathematical Biosciences and Engineering, 2022, 19(12): 11840-11853. doi: 10.3934/mbe.2022552
    [8] Wencong Zhang, Yuxi Tao, Zhanyao Huang, Yue Li, Yingjia Chen, Tengfei Song, Xiangyuan Ma, Yaqin Zhang . Multi-phase features interaction transformer network for liver tumor segmentation and microvascular invasion assessment in contrast-enhanced CT. Mathematical Biosciences and Engineering, 2024, 21(4): 5735-5761. doi: 10.3934/mbe.2024253
    [9] Yingjian Yang, Qiang Li, Yingwei Guo, Yang Liu, Xian Li, Jiaqi Guo, Wei Li, Lei Cheng, Huai Chen, Yan Kang . Lung parenchyma parameters measure of rats from pulmonary window computed tomography images based on ResU-Net model for medical respiratory researches. Mathematical Biosciences and Engineering, 2021, 18(4): 4193-4211. doi: 10.3934/mbe.2021210
    [10] Michael James Horry, Subrata Chakraborty, Biswajeet Pradhan, Maryam Fallahpoor, Hossein Chegeni, Manoranjan Paul . Factors determining generalization in deep learning models for scoring COVID-CT images. Mathematical Biosciences and Engineering, 2021, 18(6): 9264-9293. doi: 10.3934/mbe.2021456
  • We quantify a recent five-category CT histogram based classification of ground glass opacities using a dynamic mathematical model for the spatial-temporal evolution of malignant nodules. Our mathematical model takes the form of a spatially structured partial differential equation with a logistic crowding term. We present the results of extensive simulations and validate our model using patient data obtained from clinical CT images from patients with benign and malignant lesions.


    1. Introduction

    Non-small-cell lung carcinomas (NSCLC) are the most common epithelial lung cancers. The development of thin slice CT (computed tomography) scans, coupled with new recommendations for lung cancer screening in high risk patients, has led to increased detection of subtle pulmonary subsolid or nonsolid nodules in the lungs [13]. CT scan x-rays measure these nodules, also known as ground glass opacities (GGOs), as the partial filling of air spaces in the lungs by exuded fluids. Published recommendations [4], [20], [21], [22] for how to follow GGOs over time depends only on nodule size and the presence or absence of a solid component. Recent work has demonstrated the utility of volumetric CT (vCT) for diagnosis of cancer in solid nodules by measuring growth rate over time. For these cancers, which include adenocarcinoma, a growth rate given by a volume doubling time (DT) less than 400 days is predictive of malignancy [12]. However, GGOs often grow slowly in size, thus giving a high false negative rate when using nodule volume as the imaging parameter. Additionally, GGOs can be difficult to segment on CT, making assessment of growth using vCT problematic. In this work, we investigate the potential to assess GGO growth based on a quantitative change in its 3-D density histogram, irrespective of the nodule size or presence of a solid component.

    A recent report correlated five categories of CT histogram with histopathological characteristics and recurrence-free survival times [15]. Our objective is to model these five qualitative GGO measurement histogram categories, and their interpretations of tumor progression, to a quantitative dynamic mathematical model of tumor growth, which also allows estimation of tumor DT.

    The mathematical model we use for the spatial-temporal evolution of a GGO is a diffusive logistic partial differential equation. We assume cell mass grows almost exponentially in an early time phase from an initial condition consisting of a small nodule, but ultimately slows in growth as time advances. CT scans are quantified in Hounsfield units (HU), which measure radio-density. Since HU reflect tissue density as the partial filling of air spaces in the lungs by exuded fluids, it is possible for the tumor to increase in density without increasing in physical size on CT, by tumor cells gradually filling in available lung air space (see Figure 1). Therefore, visually observed CT scans may show boundaries of the tumor that do not change for a considerable amount of time, but which may increase in density. This change is reflected in the CT histogram; hence, it may be possible to quantify tumor growth based on subtle changes in the CT histogram.

    Figure 1. Photomicrograph showing a small lung area at the microscopic level. Lighter pink areas are representing the thickened alveolar walls and the darker purple ones are cancer cells lining up along the walls. As the tumor grows further, it will fill the white air spaces between the alveolar walls, thereby shifting the density histogram closer to water.

    We identify the five histogram categories formulated in [15], which are based on qualitative HU histogram signatures, with the outputs of our mathematical model, which is based on the time dependent spatial density u(t,x) of tumor cells in a spatial region Ω of the lung. The identification at a given time t is based on the fraction of values of CT scan histogram output and model output u(t,x) in specified subintervals of [0,1]. The growth and diffusion parameters in the model equation are used to identify the connection over a time series of the two outputs. The model output is then used to identify the DT values for the time series of CT scan histograms. We illustrate the usability of the model for diagnosis of lung cancer with its comparison to CT lung scan data through four clinical patient case studies.

    The Tennessee Valley Healthcare System VA Hospital Institutional Review Board approved the analysis of the anonymized CT scan data used in this paper and waived the need for informed consent.


    2. Materials and methods


    2.1. Mathematical model

    It has been recently documented that spatial intra-tumor heterogeneity plays an important role in lung cancer development at both the micro-molecular and at the macro-visible level [4], [20]. At the microscopic level and at the early stages of pulmonary adenocarcinoma in situ (previously bronchioloalveolar cell carcinoma), cancer cells align along alveolar walls in a so-called lepidic pattern. As the tumor invades the air spaces, it becomes more dense on CT.

    Mathematical models of tumor growth in spatial regions have been developed by many researchers, including [1], [3], [10], [11], [17], [18], [24], [25]. Many mathematical models have been designed specifically to connect to CT scan imaging, including [2], [5], [8], [9], [16], [26], [27]. Our goal is to develop a mathematical model that aids lung CT scan analysis, and therefore our model captures tumor spatial growth dynamics at the macro-visible level. Our model has the following form of a diffusion partial differential equation with a growth-limiting logistic term:

    ut(t,x)=(bu(t,x))+au(t,x)(1u(t,x)+ub(x)um),t>0,xΩ; (1)
    u(t,x)|Ω=0,t0;u(0,x)=u0(x),xΩ. (2)

    In the model above u(t,x) stands for the density of tumor cells at time t and spatial position xΩ, where ΩR2 is the observed physical area of the lung (typically a rectangular or disc-like area). We focus here on the 2-dimensional case, which exhibits the essential features of the underlying dynamics of lung tumor growth, and is also comparable to the clinical appearance of CT scan patient data. The 2-dimensional region ΩR2 can be viewed as representative thin slice of the tumor nodule. In future work we will consider a domain ΩR3, which is more realistic, but with numerical simulations much more time consuming. Note that for simplicity we imposed Dirichlet boundary condition in our model (1)-(2). This is because we are interested in the short term behavior of solutions, with initial tumor cell distributions supported at (or near) the center of the domain.

    The parameter a is the logistic growth rate. The parameter b is the diffusion coefficient, which determines the speed of spatial tumor propagation. In the examples the parameters a and b were qualitatively determined to match the CT histograms and the model outputs for each patient. In future work we will use formal optimization methods to specify these parameters. The initial conditions u0(x) were determined by random choice of clusters of Gaussians and then chosen for compatibility with the CT data. In future work we will develop formal methods for assigning these initial conditions specifically to patient data.

    The maximum of u(t,x) at any xΩ is umub(x). The units of u(t,x) are density units of tumor mass per unit area, which we scale to allow comparison with CT scan GGO) HU values. We take the carrying capacity parameter um=1100 as the maximum cell density at any location. We convert u(t,x) units into CT scan HU by subtracting 1000. Most body soft tissue has HU values somewhere between water (HU = 0) and blood (HU 50) due to the high iron content in blood; hence the upper limit of our histogram scale of +100 (or 1100 on the u(t, x) scale). Thus, u(t,x) values range between 0 and 1100, corresponding to HU values between 1000 and 100.

    In the subsequent section when we present our simulation results, the cell density u(t,x) will be compared to histograms represented in HU, which are volume averaged values of mixtures of air and water, HUair=1000, HUwater=0 [14]. Normal lung histograms are centered around HU=750, reflecting about 75% air. As a tumor grows, more tissue density (water) displaces the air and typically shifts the histogram to the right. Note that Hounsfield units are integer valued.

    The spatial growth of the tumor in model (1)-(2) is limited by the normal background lung cell distribution, denoted by the time-independent background density function ub(x). Tumor growth concentrates in micro-environmental regions of lung tissue vascularization [4], where the background density ub(x) is higher. The initial values and parameters are qualitatively fitted to each patient CT scan histogram data. In future work, formal optimization procedures will be developed for quantitative fittings based on large numbers of patient data.

    The proof of existence of unique solutions of model (1)-(2) is provided in the Appendix. Note that, here we are mainly interested in the early transient behavior of solutions of the model, and not the long-term asymptotic behavior.


    2.2. CT scan histogram categories

    The five CT scan histogram categories presented in [15] are summarised in Table 1 below. The classification of these categories is qualitative and subject to interpretation. The classifications of patient examples in [15] were each constructed by visual assessment of two expert observers, using a decision tree algorithm, with disagreements resolved by consensus. The histograms in the study in [15] were given in terms of continuous smoothed-out renderings of the histogram bar graphs, which allowed easier determination of category type. In our study we use actual histogram bar graphs, which preserve more information. In general, the classification of category for a given patient data set is necessarily subjective, and in fact, some patient data in our database do not readily fit any of the classifications. Our main goal is to construct a model that fits patient CT scan histogram data, rather than a model that fits the interpretation of these data according to the classification scheme in Kawata et al. We believe that our model simulations will aid in the designation of these categories for individual patients.

    Table 1. The five CT scan histogram categories.
    TypeDescription
    αhigh peak at low HU values and no peak at high HU values
    βmedium peak at low HU values and no peak at high HU values
    γlow peak at low HU values and lower peak at high HU values
    δlow peak at low HU values and higher peak at high HU values
    ϵlow peak at low HU values and very high peak at high HU values
     | Show Table
    DownLoad: CSV

    To compare model output to patient data for a time series of CT scan histograms for a given patient, we will use a quantitative determination of the fractions of both CT scan histogram outputs and model (1)-(2) outputs. The CT scan fractional histogram outputs are the fractions of histogram bar heights in a given range of HU. The fractional model outputs are the integrals of the density function u(t,x) over a given range of values, divided by the integral of all values of the density function u(t,x), with both integrals over all of Ω. We assign three output ranges as presented in Table 2 below (other choices are also possible).

    Table 2. The three output fractions.
    fractionCT scan histogram output at time tmodel output u(t,x)
    f1 <600 HU <500
    f2between 600 HU and 100 HUbetween 500 and 1000
    f3 >100 HU >1000
     | Show Table
    DownLoad: CSV

    We use the output fractions f1,f2,f3 to compare a time series of CT scan histograms for a given patient with the model (1)-(2) by specifying a time-independent background density ub(x), an initial condition u0(x) (corresponding to the baseline histogram), the logistic parameter a, and the diffusion parameter b. We then calculate the doubling time DT of the tumor from the model output.


    3. Results

    We provide here the results of simulations for four case studies, all compared to patient data. Our patient data and model simulation codes (developed in MATHEMATICA) are available upon request to the authors. All histograms, for both CT scan data and model simulations, are constructed with binning width of 10 HU per bin. We note that for each simulation the initial density u0(x) is formulated as a 2-dimensional Gaussian, and the background density ub(x) is formulated as an array of 2-dimensional Gaussians, which are parameterized so that the histogram of u0(x)+ub(x) corresponds approximately to the baseline CT histogram in each simulation. These inputs are viewed as representative of the tumor at the macro-level. In future work these inputs will be formulated at the micro-level as in Figure 1, which requires much greater detail and much more extensive computing resources for running the simulations.


    3.1. Patient 1

    Patient 1 is an example of a biopsy proven benign GGO. In Figure 2 we show CT scan images for Patient 1 at five time points. Patient 1 data consists of CT scan histograms in a series of five time points over approximately two years. These five histograms, with their category type and fractional values f1, f2, f3, are given in Figure 3. For the model simulation of Patient 1, we have taken the time points (in days) as t0=0, t1=87, t2=228, t3=643, and t4=692, corresponding to the dates in Figure 3. In Figure 4 we graph the initial tumor spatial density plus background density u0(x)+ub(x), in alignment with the CT scan histogram at baseline t=0, shown in Figure 3, and the tumor spatial density u(t,x) at t=692.

    Figure 2. Patient 1: Five serial CT images spanning 826 days (as detailed in the text) for a biopsy proven benign GGO (arrow).
    Figure 3. CT scan histograms of Patient 1. A: 9/21/12, type β, f1=0.94, f2=0.06, f3=0.0. B: 12/17/12, type β, f1=0.92, f2=0.08, f3=0.0. C: 5/7/13, type β, f1=0.94, f2=0.06, f3=0.0. D: 6/26/14, type β, f1=0.92, f2=0.08, f3=0.0. E: 8/14/14, type β, f1=0.95, f2=0.05, f3=0.0.
    Figure 4. Patient 1 model simulation. A: the initial tumor spatial density u(0,x,y). B: the initial spatial density of the tumor plus the background spatial density u(0,x,y)+ub(x,y). C: the tumor spatial density u(692,x,y) at time t=692 days.

    In Figure 5 we graph the histogram plots (with bin width 10) of the model simulation of Patient 1 at the five time points as in Figure 3, where the values of u(t,x) are shifted by 1000 to correspond to HU. The category type and fractional values f1, f2, f3 at each time point are given in the Figure 5 legend. The histograms in Figure 3 and Figure 5 show relatively good alignment, all with type β. The histogram fractions for the CT scan data and the model simulations are compared in Figure 6. From these histogram plots we see that the tumor does not progress in category type. The parameters for Patient 1 and the doubling time obtained from the simulation are shown in Table 3.

    Figure 5. Model simulation histograms of Patient 1. A: 9/21/12, type β, f1=0.97, f2=0.03, f3=0.0. B: 12/17/12, type β, f1=0.98, f2=0.02, f3=0.0. C: 5/7/13, type β, f1=0.98, f2=0.02, f3=0.0. D: 6/26/14, type β, f1=0.95, f2=0.05, f3=0.0. E: 8/14/14, type β, f1=0.94, f2=0.06, f3=0.0.
    Figure 6. Histogram fractions f1,f2,f3 of Patient 1 for CT scan data and model output.
    Table 3. Model parameters and simulation doubling times. Units of a are 1/ time units and units of b are area units2/time units.
    Patient a bDoubling time from baseline
    1 0.003 0.02353 days
    2 0.002 0.006 687 days
    3 0.004 0.001 380 days
    4 0.012 0.001 115 days
     | Show Table
    DownLoad: CSV

    3.2. Patient 2

    Patient 2 is an example of a benign GGO nodule. In Figure 7 we show CT scan images for Patient 2 at six time points. Patient 2 CT scan histograms at six time points, their category type, and fractional values f1, f2, f3, are given in Figure 8. For the model simulation of Patient 2, we have taken the six time points (in days) as t0=0, t1=107, t2=198, t3=386, t4=568, and t5=932 corresponding to the dates in Figure 8. In Figure 9 we graph the initial tumor spatial density plus background density u0(x)+ub(x), in alignment with the CT scan histogram at baseline t=0, shown in Figure 8, and the tumor spatial density u(t,x) at t=932.

    Figure 7. Patient 2: Six serial CT images over a span of 932 days for a stable GGO (arrow), clinically considered benign due to lack of change in size or density.
    Figure 8. CT scan histograms of Patient 2. A: 5/22/12, type β, f1=0.94, f2=0.06, f3=0.0. B: 9/6/12, type β, f1=0.91, f2=0.09, f3=0.0. C: 12/6/12, type β, f1=0.93, f2=0.07, f3=0.0. D: 6/12/13, type β, f1=0.92, f2=0.08, f3=0.0. E: 12/11/13, type β, f1=0.89, f2=0.11, f3=0.0. F: 12/10/14, type β, f1=0.90, f2=0.10, f3=0.0.
    Figure 9. Patient 2 model simulation: A: The initial spatial density of the tumor plus the background spatial density u(0,x,y)+ub(x,y). B: The initial tumor spatial density u(0,x,y). C: The tumor spatial density u(932,x,y) at time t=932 days.

    In Figure 10 we show the histogram plots (with bin width 10) of the model simulation of Patient 2 at the six time points as in Figure 8, where the values of u(t,x) are shifted again by 1000 so that they correspond to HU. The category type and fractional values f1, f2, f3 at each time point are given in the Figure 10 legend. The histograms in Figure 8 and Figure 10 show relatively good alignment, all with type β. The histogram fractions for the CT scan data and the model simulations are compared in Figure 11. From these histogram plots we see that the tumor does not progress in category type. The parameters for Patient 2 and the doubling times obtained from the simulation are shown in Table 3. The doubling time obtained from the simulation is in the range of a benign nodule.

    Figure 10. Model simulation histograms of Patient 2. A: 5/22/12, type β, f1=0.92, f2=0.08, f3=0.0. B: 9/6/12, type β, f1=0.91, f2=0.09, f3=0.0. C: 12/6/12, type β, f1=0.91, f2=0.09, f3=0.0. D: 6/12/13, type β, f1=0.89, f2=0.11, f3=0.0. E: 12/11/13, type β, f1=0.87, f2=0.13, f3=0.0. F: 12/10/14, type β, f1=0.84, f2=0.16, f3=0.0.
    Figure 11. Histogram fractions f1,f2,f3 of Patient 2 for CT scan data and model output.

    3.3. Patient 3

    Patient 3 is an example of atypical cells highly suspicious for adenocarcinoma by biopsy. In Figure 12 we show CT scan images for Patient 3 at four time points. Patient 3 CT scan histograms (with bin width 10 HU) at the four time points, with their category type and fractional values f1, f2, f3 are shown in Figure 13. For the model simulation of Patient 3, we have taken the time points (in days) as t0=0, t1=574, t2=826, t3=917 (corresponding to the dates in Figure 13), and two additional time points beyond the data times as t4=1,217 and t5=1,517. In Figure 14 we graph the initial tumor spatial density plus background density u0(x)+ub(x), in alignment with the CT scan histogram at baseline t=0, shown in Figure 13, and the tumor spatial density u(t,x) at t=917.

    Figure 12. Patient 3: Four serial CT images spanning 917 days for atypical cells (arrow) highly suspicious for adenocarcinoma by biopsy.
    Figure 13. CT scan histograms of Patient 3. A: 10/20/10, type β, f1=0.74, f2=0.23, f3=0.03. B: 5/16/11, type β, f1=0.69, f2=0.24, f3=0.07. C: 1/23/13, type γ, f1=0.69, f2=0.22, f3=0.09. D: 4/24/13, type γ, f1=0.63, f2=0.24, f3=0.13.
    Figure 14. Patient 3 model simulation: A: The initial spatial density of the tumor plus the background spatial density u(0,x,y)+ub(x,y). B: The initial tumor spatial density u(0,x,y). C: The tumor spatial density u(917,x,y) at time t=917 days.

    In Figure 15 we show the histogram plots of the model simulation for Patient 3 at the six time points as shown in Figure 13, where the values of u(t,x) are shifted by 1000 to correspond to HU. The category type and fractional values f1, f2, f3 at each time point are given in the Figure 15 legend. The first four histograms in Figure 13 and Figure 15 show relatively good alignment, with type progression from β to γ. Through the two additional time points in the simulation we see the progression of the tumor through type γ. The histogram fractions for the CT scan data and the model simulations are compared in Figure 16. The parameters for Patient 3 and the doubling time obtained from the simulation are shown in Table 3. The doubling time obtained from the simulation is in the range of non-small cell lung cancer.

    Figure 15. Model simulation histograms of Patient 3. A: 10/20/10, type β, f1=0.78, f2=0.22, f3=0.0. B: 5/16/12, type β, f1=0.62, f2=0.38, f3=0.0. C: 1/23/13, type γ, f1=0.55, f2=0.42, f3=0.03. D: 4/24/13, type γ, f1=0.52, f2=0.43, f3=0.05. E: 4/24/13 + 300 days, type δ, f1=0.44, f2=0.43, f3=0.13. F: 4/24/13 + 600 days, type δ, f1=0.37, f2=0.40, f3=0.23.
    Figure 16. Histogram fractions f1,f2,f3 of Patient 3 for CT scan data and model output.

    3.4. Patient 4

    Patient 4 is an example of a proven adenocarcinoma that started as a GGO that increased in density on CT over time. In Figure 17 we show CT scan images for Patient 4 at four time points. Patient 4 CT scan histograms (with bin width 10 HU) at four time points, with their category type and fractional values f1, f2, f3 are shown in Figure 18. For the model simulation of Patient 4, we have taken the time points (in days) as t0=0, t1=239, t2=423, t3=471, and two additional time points beyond the data times as t4=600 and t5=750. In Figure 19 we show the graph of the initial spatial density u0(x), the background spatial density ub(x), and their sum, in alignment with the CT scan histogram at baseline t=0 shown in Figure 18.

    Figure 17. Patient 4: Four CT image recordings of a suspicious nodule spanning 471 days (arrow).
    Figure 18. CT scan histograms of Patient 4. A: 10/2/13, type β, f1=0.72, f2=0.24, f3=0.04. B: 5/28/14, type γ, f1=0.54, f2=0.35, f3=0.11. C: 11/28/14, type γ, f1=0.56, f2=0.33, f3=0.11. D:1/15/15, type γ, f1=0.46, f2=0.36, f3=0.18.
    Figure 19. Patient 4 model simulation: A: The initial spatial density of the tumor plus the background spatial density u(0,x,y)+ub(x,y). B: The initial tumor spatial density u(0,x,y). C: The tumor spatial density u(471,x,y) at time t=471 days.

    In Figure 20 we show the histogram plots of the model simulation for Patient 4 at the four time points as shown in Figure 18, where the values of u(t,x) are shifted by 1000 to correspond to HU. The category type and fractional values f1, f2, f3 at each time point are given in the Figure 20 legend. The first four histograms in Figure 17 and Figure 20 show relatively good alignment, with type progression from β to γ. The histogram fractions for the CT scan data and the model simulations are compared in Figure 21. Through the two additional time points we see the progression of the tumor through type δ to type ϵ. The parameters for Patient 4 and the doubling time obtained from the simulation are shown in Table 3. The doubling time obtained from the simulation is in the range of non-small cell lung cancer.

    Figure 20. Model simulation histograms of Patient 4. A: 9/21/12, type β, f1=.90, f2=0.10, f3=0.0. B: 12/17/12, type γ, f1=0.66, f2=0.34, f3=0.0. C: 5/7/13, type γ, f1=0.47, f2=0.40, f3=0.13. D: 6/26/14, type γ, f1=0.43, f2=0.39, f3=0.17. E: 8/14/14, type δ, f1=0.33, f2=0.37, f3=0.30. F: 10/20/15, type ϵ, f1=0.23, f2=0.33, f3=0.44.
    Figure 21. Histogram fractions f1,f2,f3 of Patient 4 for CT scan data and model output.

    In Figure 22 we graph the total tumor mass from the model simulations for each patient over time, where mass is scaled to 1.0 at time 0. Patients 1 and 2 have smaller growth than Patients 3 and 4, corresponding to their smaller growth parameter a. The model simulations allow calculation of the tumor doubling times, as well as tracking of the tumor growth over the span of CT scan time series. The model simulations also allow growth projections for additional times beyond the CT scan data, as demonstrated for Patients 3 and 4 (see also Figure 22).

    Figure 22. Total tumor mass growth curves from model simulations. Black dots are time points corresponding to CT scan data for patients 1, 2, 3, 4. Red dots are for two additional time points for Patients 3 and 4. The values are scaled to 1.0 at time 0.

    4. Discussion

    In a recent paper [15], a qualitative five-category classification method was proposed for analyzing NSCLC, and its utility justified using statistical tools. The results indicated a satisfactory inter-observer agreement simply through visual assessment of CT histograms. Our goal here has been to quantify the five categories in [15] in terms of a dynamic spatial model of tumor growth; and to connect the temporal dynamics of the categories to tumor DT. We have compared CT scan data and model outputs for four patient studies. For each patient, we see good agreement between these data and model outputs, in terms histogram categories and HU fractional ranges.

    In the current work we hypothesized that the five categories identified in [15] actually correspond to temporal tumor progression. Indeed, Kawata [15] already speculated that change from type α to β and from β to γ may indicate tumor progression.

    Our results show that model (1)-(2) supports the five category classification in adenocarcinoma in situ. Further, these five categories can be viewed as a hypothesized 5-step lung cancer progression theory. Moreover, since it takes into account the spatial heterogeneity of the tumor, which is particularly important for irregular nodules investigated here, the model gives us a tool to estimate tumor mass doubling times using CT histogram data only.

    Major challenges for application of the model (1)-(2)are the identifications of the initial tumor nodule characteristics, the background non-tumor bias parameter ub(x), the carrying capacity parameter um, the spatial diffusion parameter b, and tumor growth parameter a. Our goal here has been to demonstrate that model (1)-(2) does correlate well with tumor growth data given by CT scan data represented with GGO histograms. Formal procedures to quantify these identifications of initial data and parameters for general patient data will be carried out in future work.


    5. Outlook

    Our model already shows very good agreement with patient data, and the 5-category classification of GGOs. Future improvements of the mathematical model may involve:

    ● Full 3-dimensional simulations.

    ● Systematically analyze the simulation outcomes as functions of the model parameters and initial condition (transient vs asymptotic behavior, is there a globally stable steady state?).

    ● Inclusion of nonlinear diffusion to account for a more realistic description of tumor spatial growth (in particular to model competition effects).

    ● To include different type of placement processes for the tumor cells (other than diffusion) to account for the complex spatial structure of the lung.

    Estimation of tumor doubling time in GGOs has not been described. This work offers a method to compute growth rate of GGOs as a predictive biomarker of malignancy, similar to that used for solid nodules using volumetric CT. Further work is needed to investigate the impact of different reconstruction algorithms and reconstructed image quality on the estimate of GGO growth rate.


    6. Supporting information

    Global behaviour of solutions. The basic mathematical theory of general classes of nonlinear reaction diffusion equations of the type (1)-(2) is well understood. However, for completeness, here we provide a concise proof of the global existence and positivity of solutions of our model in the biologically relevant state space of Lebesgue integrable functions L1(Ω)=:X. In particular, to establish the global existence of mild solutions we implement a framework as in [23] for a structured population model, see also [19].

    We set K:=L1+(Ω) (the positive cone of L1, which is closed) and we recast model (1)-(2) in the form of a semilinear abstract Cauchy problem as follows.

    dudt=Au+F(u),u(0)=u0K, (3)

    where

    Au=(bu)+au(1ubum), (4)
    D(A)={vW2,1(Ω)|v(x)=0xΩ}, (5)
    F(u)=au2um. (6)

    We say that the abstract semilinear problem (3) satisfies the sub-tangential condition (see e.g. [23]) with respect to K, if

    limh0+d(K,T(h)u+hF(u))h=0, (7)

    where T is the linear semigroup generated by A, and d is the usual distance function. We also recall the notation (,) introduced for a semi-inner product on X. Below X denotes the dual of the Banach space X, and (,) the natural pairing between elements of X and X.

    (u,v):=minvX{(u,v)|||v||=||v||,(v,v)=||v||2}.

    We recall the following result from [23], see also [19].

    Let X,K,A and F as defined above, and assume that F is locally Lipschitz and bounded. Further assume that the sub-tangential condition (7) holds, and that there exist ω,κR such that (Au,u)ω|u|2, for all uD(A); and (F(u),u)κ|u|2, for all uK. Then, for each u0K, there exists a unique mild solution u(t) to (3) for all t>0.

    We now apply this result to our model (1)-(2).

    Theorem 6.1. Assume that ubC1(¯Ω), and a,b>0. Then, for any initial condition u0K, model (1)-(2) admits a mild (semigroup) solution u(t)K, for all times t>0.

    Proof. It follows from the assumptions that the densely defined operator A defined in (4)-(5) generates a positive strongly continuous semigroup T(t) on L1(Ω). Note that the nonlinear operator F cannot be defined on the whole state space X, but F is locally Lipschitz and maps bounded sets BK into bounded sets F(B). To establish the global existence of solutions, note that in our situation since T leaves K invariant, the sub-tangential condition (7) simplifies as follows (see also Lemma C in [23]).

    limh0+d(K,u+hF(u))h=0, (8)

    which is easily seen to hold true, as for all uK we have F(u)<0.

    Next note that in our setting we have

    (F(u),u)=minuL(Ω){aumΩu2u|||u||1=||u||,Ωuu=(Ω|u|)2}. (9)

    Hence for every uK we may take u||u||1=Ωu (constant function), which shows that (F(u),u)0 holds. Finally, note that (Au,u)ω|u|2, for all uD(A) holds with ω:=s(A)<, the spectral bound of A.

    Our model (1)-(2) always admits the trivial steady state u0. For a large enough, the existence of a strictly positive steady state can be established using the general framework developed in [6]. In particular, we can define a parametrized family of linear operators as follows:

    Φvu=(bu)+au(1ubum)auumv, (10)
    D(Φv)={uW2,1(Ω)|u(x)=0xΩ},vK. (11)

    It is then shown that for a large enough, s(Φ0)>0 holds, and that the function defined as f:αs(Φαv) is monotone decreasing for every vK. This then allows one to define a fixed point map on the level set S:={vK|s(Φv)=0}, which yields the existence of a positive steady state of (1)-(2), see [6] for more details.

    We also note that applying earlier results by Cantrell and Cosner from [7] (see in particular Theorem 3.1 in [7]) would also allow us to obtain sufficient conditions for the existence of a globally stable unique positive steady state for a large enough.


    Acknowledgments

    We acknowledge support provided by the National Cancer Institute U01CA152662. We also thank the Royal Society, which supported visits of József Farkas to Vanderbilt University.


    [1] [ D. Ambrosi,F. Mollica, On the mechanics of a growing tumor, Int. J. Eng. Sci., 40 (2002): 1297-1316.
    [2] [ H. Ammari, Mathematical Modeling in Biomedical Imaging 1. Electrical and Ultrasound Tomographies, Anomaly Detection, and Brain Imaging, Springer Science and Business Media, New York, 2009.
    [3] [ A. R. A. Anderson,A. M. Weaver,P. T. Cummings,V. Quaranta, Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment, Cell, 127 (2006): 905-915.
    [4] [ F. R. Balkwill,M. Capasso,T. Hagemann, The tumor microenvironment at a glance, J. Cell Sci., 125 (2012): 5591-5596.
    [5] [ T. M. Buzug, Computed Tomography, from Photon Statistics to Modern Cone-beam CT, Springer-Verlag, Berlin-Heidelberg-New York, 2008.
    [6] [ Á. Calsina,J. Z. Farkas, Positive steady states of nonlinear evolution equations with finite dimensional nonlinearities, SIAM J. Math. Anal., 46 (2014): 1406-1426.
    [7] [ R. S. Cantrell,C. Cosner, Diffusive logistic equations with indefinite weights: Population models in disrupted environments Ⅱ, SIAM J. Math. Anal., 22 (1991): 1043-1064.
    [8] [ O. Clatz,M. Sermesant,P.-Y. Bondiau,H. Delingette,S. K. Warfield,G. Malandain,N. Ayache, Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation, IEEE Trans. Med. Imag., 24 (2005): 1334-1346.
    [9] [ F. Cornelis,O. Saut,P. Cumsille,D. Lombardi,A. Iollo,J. Palussiere,T. Colin, In vivo mathematical modeling of tumor growth from imaging data: Soon to come in the future?, Diag. Inter. Imag., 94 (2013): 593-600.
    [10] [ H. Enderling,M. A. J. Chaplain, Mathematical modeling of tumor growth and treatment, Curr. Pharm. Des., 20 (2014): 4934-4940.
    [11] [ R. A. Gatenby,P. K. Maini,E. T. Gawlinski, Analysis of a tumor as an inverse problem provides a novel theoretical framework for understanding tumor biology and therapy, Appl. Math. Lett., 15 (2002): 339-345.
    [12] [ C. I. Henschke,D. F. Yankelevitz,R. Yip,A. P. Reeves,D. Xu,J. P. Smith,D. M. Libby,M. W. Pasmantier,O. S. Miettinen, Lung cancers diagnosed at annual CT screening: Volume doubling times, Radiology, 263 (2012): 578-583.
    [13] [ C. I. Henschke,R. Yip,J. P. Smith,A. S. Wolf,R. M. Flores,M. Liang,M. M. Salvatore,Y. Liu,D. M. Xu,D. F. Yankelevitz, CT screening for lung cancer: Part-solid nodules in baseline and annual repeat rounds, Am. J. Roentgenol, 207 (2016): 1176-1184.
    [14] [ G. N. Hounsfield, Computed medical imaging, Nobel Lecture, J. Comput. Assist. Tomogr., 4 (1980): 665-674.
    [15] [ Y. Kawata,N. Niki,H. Ohmatsu,M. Kusumoto,T. Tsuchida,K. Eguchi,M. Kaneko,N. Moriyama, Quantitative classification based on CT histogram analysis of non-small cell lung cancer: Correlation with histopathological characteristics and recurrence-free survival, Med. Phys., 39 (2012): 988-1000.
    [16] [ E. Konukoglu,O. Clatz,B. H. Menze,B. Stieltjes,M-A. Weber,E. Mandonnet,H. Delingette,N. Ayache, Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations, IEEE Trans. Med. Imag., 29 (2010): 77-95.
    [17] [ Y. Kuang, J. D. Nagy and S. E. Eikenberry, Introduction to Mathematical Oncology, Mathematical and Computational Biology Series, Taylor & Francis Group, Boca Raton-London-New York, 2016.
    [18] [ J. S. Lowengrub,H. B. Feiboes,F. Jin,Y.-I. Chuang,X. Li,P. Macklin,S. M. Wise,V. Cristini, Nonlinear modelling of cancer: Bridging the gap between cells and tumors, Nonlinearity, 23 (2010): R1-R91.
    [19] [ R. H. Martin, Nonlinear Operators and Differential Equations in Banach Spaces, Pure and Applied Mathematics, Wiley-Interscience, John Wiley & Sons, New York-London-Sydney, 1976.
    [20] [ D. Morgensztern,K. Politi,R. S. Herbst, EGFR Mutations in non-small-cell lung cancer: Find, divide, and conquer, JAMA Oncol., 1 (2015): 146-148.
    [21] [ D. P. Naidich,A. A. Bankier,H. MacMahon,C. M. Schaefer-Prokop,M. Pistolesi,J. M. Goo,P. Macchiarini,J. D. Crapo,C. J. Herold,J. H. Austin,W. D. Travis, Recommendations for the management of subsolid pulmonary nodules detected at CT: A statement from the Fleischner Society, Radiology, 266 (2013): 304-317.
    [22] [ National lung screening trial research team, Reduced lung-cancer mortality with low-dose computed tomographic screening, N. Engl. J. Med., 365 (2011), 395–409.
    [23] [ J. Prüss, Equilibrium solutions of age-specific population dynamics of several species, J. Math. Biol., 11 (1981): 65-84.
    [24] [ R. Rockne, E. C. Alvord, Jr., M. Szeto, S. Gu, G. Chakraborty and K. R. Swanson, Modeling diffusely invading brain tumors: An individualized approach to quantifying glioma evolution and response to therapy, in Selected Topics in Cancer Modeling: Genesis, Evolution, Immune Competition and Therapy, Modeling and Simulation in Science, Engineering and Technology Series, Birkh¨auserBoston, Boston, MA, 2008,207–221.
    [25] [ K. R. Swanson,C. Bridge,J. D. Murray,E. C. Alvord Jr., Virtual and real brain tumors: Using mathematical modeling to quantify glioma growth and invasion, J. Neurol. Sci., 216 (2003): 1-10.
    [26] [ C. H. Wang,J. K. Rockhill,M. Mrugala,D. L. Peacock,A. Lai,K. Jusenius,J. M. Wardlaw,T. Cloughesy,A. M. Spence,R. Rockne,E. C. Alvord Jr.,K. R. Swanson, Prognostic significance of growth kinetics in newly diagnosed glioblastomas revealed by combining serial imaging with a novel biomathematical model, Can. Res., 69 (2009): 9133-9140.
    [27] [ A. Y. Yakovlev, A. V. Zorin and B. I. Grudinko, Computer Simulation in Cell Radiobiology, Lecture Notes in Biomathematics, 74, Springer-Verlag, Berlin-Heidelberg-New York, 1988.
  • This article has been cited by:

    1. Qi Mao, Shuguang Zhao, Lijia Ren, Zhiwei Li, Dongbing Tong, Xing Yuan, Haibo Li, Intelligent immune clonal optimization algorithm for pulmonary nodule classification, 2021, 18, 1551-0018, 4146, 10.3934/mbe.2021208
    2. Adriana Del Pino Herrera, Meghan C. Ferrall-Fairbanks, A war on many fronts: cross disciplinary approaches for novel cancer treatment strategies, 2024, 15, 1664-8021, 10.3389/fgene.2024.1383676
  • Reader Comments
  • © 2018 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4228) PDF downloads(680) Cited by(2)

Article outline

Figures and Tables

Figures(22)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog