Mathematical modeling of multiple pathways in colorectal carcinogenesis using dynamical systems with Kronecker structure

Abstract

Like many other types of cancer, colorectal cancer (CRC) develops through multiple pathways of carcinogenesis. This is also true for colorectal carcinogenesis in Lynch syndrome (LS), the most common inherited CRC syndrome. However, a comprehensive understanding of the distribution of these pathways of carcinogenesis, which allows for tailored clinical treatment and even prevention, is still lacking. We suggest a linear dynamical system modeling the evolution of different pathways of colorectal carcinogenesis based on the involved driver mutations. The model consists of different components accounting for independent and dependent mutational processes. We define the driver gene mutation graphs and combine them using the Cartesian graph product. This leads to matrix components built by the Kronecker sum and product of the adjacency matrices of the gene mutation graphs enabling a thorough mathematical analysis and medical interpretation. Using the Kronecker structure, we developed a mathematical model which we applied exemplarily to the three pathways of colorectal carcinogenesis in LS. Beside a pathogenic germline variant in one of the DNA mismatch repair (MMR) genes, driver mutations in APC, CTNNB1, KRAS and TP53 are considered. We exemplarily incorporate mutational dependencies, such as increased point mutation rates after MMR deficiency, and based on recent experimental data, biallelic somatic CTNNB1 mutations as common drivers of LS-associated CRCs. With the model and parameter choice, we obtained simulation results that are in concordance with clinical observations. These include the evolution of MMR-deficient crypts as early precursors in LS carcinogenesis and the influence of variants in MMR genes thereon. The proportions of MMR-deficient and MMR-proficient APC-inactivated crypts as first measure for the distribution among the pathways in LS-associated colorectal carcinogenesis are compatible with clinical observations. The approach provides a modular framework for modeling multiple pathways of carcinogenesis yielding promising results in concordance with clinical observations in LS CRCs.

Fig 1. From the medical hypothesis over the modeling approach to the mathematical structure.

The medical hypothesis of multiple pathways of carcinogenesis is widely known for various types of cancer. (A) We present a model for this phenomenon at the example of Lynch syndrome, the most common inherited CRC syndrome, with specific key driver events in the MMR genes, CTNNB1, APC, KRAS and TP53. (B) This current medical understanding of carcinogenesis is translated into a mathematical model using a specific dynamical system, which can be represented by a graph structure, where each vertex in the graph represents a genotypic state and the edges correspond to the transition probabilities between those states. Starting with all colonic crypts in the state of all genes being wild-type and a single MMR germline variant due to Lynch syndrome, we are interested in the distribution of the crypts among the graph at different ages of the patient in order to obtain estimates for the number of crypts in specific states, e.g., adenomatous or cancerous states. (C) The underlying matrix of the dynamical system makes use of the Kronecker sum and product. It is a sparse upper triangular matrix accounting for the assumption that mutations cannot be reverted. This allows fast numerical solving by using the matrix exponential. Each nonzero entry of the matrix represents a connection between genotypic states in the graph.

Fig 2. Gene mutation graphs for independent mutational processes.

These graphs represent the possible mutation status, i.e., which mutations the alleles of the gene can have accumulated, as vertices ∅, m, l, mm, ll and ml. They are given for (A) MMR gene mutations, (B) APC mutations, (C) KRAS mutations, (D) CTNNB1 mutations, and (E) TP53 mutations. The edges connecting different vertices represent mutations, whereas self-loops, i.e., edges that connect a vertex with itself, describe no mutation occurring at the current point in time. The edges are labeled by the amount of change which happens at each point in time. Note that in the colon, biallelic mutations of CTNNB1 seem to be required to mediate an oncogenic driver effect [61, 62], leading to a gene mutation graph similar to that of APC and TP53.

Fig 3. Sparse matrix structure.

(A) The system matrix (A + B + C + D + E + F) of the linear model is a very sparse matrix, i.e., only a few entries are nonzero. These nonzero entries are colored red in the plot, which also illustrates the fact that (A + B + C + D + E + F) is an upper triangular matrix. (B) The sparsity structure of the matrix expm(A + B + C + D + E + F), which is reminiscent of a Sierpiński fractal, is due to the individual matrices being the Kronecker product and sum of matrices. The two plots also illustrate nicely how modeling sparse local interactions in the matrix (A + B + C + D + E + F) can have a more global effect in expm(A + B + C + D + E + F).

Fig 4. Gene mutation graph of APC for increasing the point mutation rate of APC after MMR deficiency.

Fig 5. Model component for the positive association of MLH1 and CTNNB1.

Part of the combined gene mutation graph for CTNNB1 and MLH1 of the matrix C. The gene mutation graphs for the other possible gene states MLH1 ∈ {l, ll}, CTNNB1 ∈ {m, ml} are defined in an analogous way.

Fig 6. Model component for increasing the LOH rate of MMR, CTNNB1 and TP53 by a factor δ + 1 after APC inactivation.

Gene mutation graph for both genes, CTNNB1 and TP53, of the component D. The gene mutation graph for MMR is defined in an analogous way.

Fig 7. Model component for the mutual enhancement of two dependencies by a factor δr_effLOH.

Part of the gene mutation graph for CTNNB1 and MLH1 after APC inactivation considered by the component E. The gene mutation graphs for the other possible gene states MLH1 ∈ {l, ll}, CTNNB1 ∈ {m, ml}, APC ∈ {ml} are defined in an analogous way.

Fig 8. Model component for increasing the mutation rate of KRAS after MMR deficiency.

Gene mutation graph of KRAS for the matrix F with the KRAS mutation rate increased by a factor ζ.

Fig 9. Number of MMR-deficient crypts over the life of a typical Lynch syndrome patient for MLH1 and MSH2.

The parameters in the model are set in such a way that the simulation results are in concordance with published data [80]. In our model, differences among genes are due to differences in coding region and gene lengths as well as the magnitude of the effects of the dependent mutational processes.

Fig 10. Number of crypts over time in a typical MLH1 carrier in combined states, like early adenomatous, advanced adenomatous and cancerous states as defined in the text for the given parameter set.

Due to the model components accounting for different genetic dependencies, the distribution of MMR-deficient and MMR-proficient, as well as the contribution of APC and CTNNB1 change for the different states. Due to the lack of suitable medical data, parameter learning was not performed in a rigorous way. As soon as data are available, this can be done using different mathematical techniques.

Fig 11. Proportion of MMR-proficient and MMR-deficient crypts in a typical MLH1 carrier in different states corresponding to the states in the classical adenoma-carcinoma sequence by Vogelstein [7].

Among the APC-/- crypts (left), the number of MMR-deficient crypts is up to 20% higher than the number of MMR-proficient ones. This difference largely increases with the subsequent KRAS activation (KRAS+) (middle) and TP53 inactivation (TP53-/-) (right) leading to the fact that almost all crypts in the last state, corresponding to a cancerous state, are MMR-deficient. These simulation results are in concordance with available data with a slight underestimation of MMR-deficient APC-/- crypts [22].

Fig 12. Comparison of MMR-deficient crypts in Lynch-like and Lynch syndrome individuals.

The number of MMR-deficient crypts is significantly higher in Lynch syndrome individuals compared to Lynch-like individuals, which matches the findings in [80].

Fig 13. Comparison of APC-/- crypts in the sporadic case and in FAP individuals, where we changed the initial value of the dynamical system as well as n_hs(APC) = 600 for FAP.

Our simulation results yield numbers below estimates found in the literature [–80]. With improved measurements, future work will adapt the parameters accordingly.