A physical mechanism of cancer heterogeneity

We studied a core cancer gene regulatory network motif to uncover possible source of cancer heterogeneity from epigenetic sources. When the time scale of the protein regulation to the gene is faster compared to the protein synthesis and degradation (adiabatic regime), normal state, cancer state and an intermediate premalignant state emerge. Due to the epigenetics such as DNA methylation and histone remodification, the time scale of the protein regulation to the gene can be slower or comparable to the protein synthesis and degradation (non-adiabatic regime). In this case, many more states emerge as possible phenotype alternations. This gives the origin of the heterogeneity. The cancer heterogeneity is reflected from the emergence of more phenotypic states, larger protein concentration fluctuations, wider kinetic distributions and multiplicity of kinetic paths from normal to cancer state, higher energy cost per gene switching, and weaker stability.


Details of the network motif
For the self activating and mutually repressing network motif, the two genes Gene A and Gene B each has two binding sites. The first binding site of gene A(B) can be bound to a monomer produced by gene B(A) and the synthesis rate of protein A(B) will be repressed by a factor λ R . The second binding site of gene A(B) can be bound to a tetramer produced by gene A (B) and the synthesis rate of A(B) will be raised by a factor λ A . We are using multiplicative model with activation and repression effectively measured as multiplicative factors instead of conventional additive factors. We believe multiplicative model is more close in practice as to how gene regulatory network makes biological logic decisions. The synthesis rates for different bound/unbound states are g 00 , g 10 = λ R g 00 , g 01 = λ A g 00 , g 11 = λ R λ A g 00 . The degradation rate for both protein A and protein B is set to be k A = k B = k. The unbinding rate for all the 4 binding sites are set equat GeneA 0α + (n + 1)B GeneB α1 + nB nA k g (n + 1)A nB k g (n + 1)B In above reactions, α represents an arbitrary state of the binding site, it can be 0 or 1. At large volume limit, protein concentration x = n V becomes continuous variable. For simplicity we absorb the volume 'V' into g, f and Xeqs. We set synthesis rates g 00 = 5, λ A = 8, λ R = 0.2, degradation rate k = 0.1, unbinding rate f = k · ω and binding rate h 1A = f P is a 16 component vector whose component P s (x, t) represents the probability of the system being at gene state s with protein concentration x at time t.
H 0 describes protein synthesis and degradation processes. It is diagonal with 16 diagonal elements. Each describes a continuous landscape at corresponding discrete gene state.
Each diagonal element L ijkl (each index i, j, k, l can be either 0 or 1) is a 'Focker-Planck' operator: On the other hand, H b describes binding/unbinding processes and is non-diagonal. It describes the 'coupling' between discrete gene states.
The diagonal elements r ijkl in H B are escaping rates from gene state ijkl.
Adiabatic Limit from Analytical Approach The gene network dynamics can be simulated with the Gillespie algorithm. 1 Various analytical approaches can be developed to approximate the adiabatic fast regulation and nonadiabatic slow regulation limit. When ω is large, we are at adiabatic limit where adiabatic approximation of fast regulation of protein to the gene compared to the protein synthesis/degradation is valid. 2,3 When H b terms representing binding/unbinding or gene switches dominate over H 0 terms representing protein synthesis/degradation, the dynamics of the fast degree of freedom (gene switching) quickly reaches equilibrium which is the steady state of H b : here ξ is a 16 components state vector and is normalized as 1 T · ξ = i ξ i = 1. Being close to ξ, the total probability distribution can be written as: where ρ is a scaler function of x and t, is a small deviation from ρξ.
The probability of the gene network with x protein concentration is the sum over the probabilities of the network at x protein concentration but being at different gene states .
That is 1 T · P = i P i ≈ ρ. Multiply the master equation using 1 T from the left, and notice we always have 1 T H b = 0. The master equation becomes: The right hand side is the sum of 16 Fokker-Planck equations result in the final form of 2-dimensional Fokker-Planck equation. The driving force has the form of Hill functions and intrinsic fluctuation is x dependent.
with driving force and diffusion coefficients: When the whole system reaches steady state (protein concentrations in addition to the gene switchings), we can quantify the landscape as : The which is similar to constant diffusion case:. The driving force can be decomposed into a gradient part plus a curl part: Non-zero flux is a measure of how far the system deviates from the detailed balance (how far away the system is from equilibrium). The non-zero flux is the origin of the irreversibility and entropy production. 4,5 As shown in Fig. 2 Path integral tells us that, among all the possible transition paths connecting initial (normal) state and final (cancer) state, the possibility of a single path is proportional to the exponential of the negative action S = The optimal path is the one that minimizes action S. 8 As has been pointed out, by using kinetic equation and energy conservation along optimal path, the action can be simplified into a line integral along the path 9,10 F l is the force projected to the path line modularized by the diffusion. The second term in the action suggests the path is irreversible since the force is not a pure gradient. The curl flux component of the force will lead to non-trivial contribution to the dynamics. This confirms the irreversibility of cancer model: the reversed optimal path will no longer be the optimal path as the new action is larger and the probability for the reversed path is smaller.

Non-Adiabatic Limit from Analytical Approach
As ω decreases, binding/unbinding processes are less frequent and we reach the non-adiabatic regime. Adiabatic approximation is no longer valid. We have to use 16 dimensional master equation explicitly.
In the small ω, or weak coupling limit. H b coupling terms are much smaller than H 0 terms. System tends to stay in one of the 16 states and rarely jumps to another. In this limit, besides normal state and cancer state, we are able to identify more intermediate states.
Especially the 'off-off' state where both genes are repressed. These states have a great impact on the cancer transition behaviour. There are also smaller peaks emerging around the cancer state, suggesting possible phenotype alternations even after the system reaches cancer state.
The moderate ω non-adiabatic region is of special interest. Conventionally in this regima solid theoretical framework and analytical results are hard to achieve. Recent progress 6,[11][12][13] showed that there's similarity between coupled Fokker-Planck equation with driving force and diffusion coefficients: The steady state probability distribution provides a global landscape picture. Compared In particular, the optimal transition path can be quantified under this unified framework in continuous extended space. Fig. 5(b) shows the optimal path obtained using minimal action approach under 6-dimensional unified landscape projected to 2-dimensional n A − n B space. Compared with optimal path at adiabatic limit, it's more tilted towards the 'off-off' state. This is due to the location change of intermediate state as a result of the weakened coupling, the change of gradient force of the landscape and non-equilibrium flux.

Stability and Heterogeinity
As we can see from Fig.2 in main text, the landscape at adiabatic limit, extreme nonadiabatic limit and moderate non-adiabatic regime are very different. In the extreme nonadiabatic region, H b terms are much smaller than H 0 terms. The gene states defined by H 0 are decoupled and we can expect as many discrete states as the dimension of H 0 matrix.
In our examples, it is up to 16 states. In the other limit where H b dominates over H 0 , adiabatic approximation is valid and the system is mapped to an N-dimensional network whose number of stable states tends to be less. The moderate non-adiabatic regime, as we can see, is different from either limit. The location of the stable states are close to the ones in the extreme non-adiabatic regime, but the smaller basins that are close to larger basins merge to major peaks which is similar to the case in adiabatic limit.
To measure the heterogeneity, we take a look at Fano factor. Fano factor is defined as variance/mean and it measures the degrees of fluctuations. It is one if the process is purely random and the distribution is Poisson. Large Fano factor indicates large deviation from Poisson distribution and large fluctuation. As Fig 3 shows, Fano factor increases as adiabaticity decreases when ω is large and doesn't change much when ω is small. This is to be expected. When ω is large the hopping processes reach equilibrium, the 16 discrete states are strongly coupled. We are left with less steady states. Fano factor decreases as ω increases in this regime. However when ω is small enough and all the discrete state peaks are decoupled, the topology of the landscape won't change and Fano factor won't change significantly. It will even decrease a little bit as ω decreases due to less fluctuations in hopping processes.
The fact that at non-adiabatic regime there are more steady states and the steady states are more fluctuating all suggest there are more phenotype states at non-adiabatic regime.
Mean First Passage Time (MFPT) is the average of first passage time. It measures the average time a successful transition takes and the transition rate is often defined as r = 1/M F P T . Fig. 7 shows how MFPT from normal state to cancer state changes with respect to adiabaticity. It suggests transitions at adiabatic limit and extreme non-adiabatic limit are rare and there exists an optimal transition rate at moderate ω non-adiabatic region.
MFPT becomes larger at small ω non-adiabatic region, since the 16 discrete landscapes are decoupled and the system tends to stay in one of the landscape for a long time. It takes several jumps from the initial state to the final state. In our example from the initial 0110 state to the final 1001 state (1 stands for bound state and 0 unbound state for each of the four binding sites), it takes at least 4 jumps and the chance of this happening is relatively rare. At adiabatic limit the hopping processes occur so frequently that the genes won't stay in one of the discrete states for a long time. There is no sufficient time for protein copy numbers to reach the transcription level corresponding to the gene state. The moderate ω non-adiabatic regime, benefits from both proper jumping probability as well as sufficient residence time for transcription and as a result has the larger transition rate. Fig. 5 in main text shows the optimal path under adiabatic limit using adiabatic approximation and moderate ω non-adiabatic regime using unified landscape approximation. The two optimal paths deviate from each other. There's no doubt that the optimal path is the transition path with maximum possibility among all the transition paths, but will optimal path give major contribution to total transition rate? Probably not in all circumstances.
We can answer this question by checking the distribution of FPT. As Fig. 4  All these results suggest moderate ω non-adiabatic regime to have larger heterogeneity.
But where does the heterogeneity come from? The number of phenotype states decreasing as ω increasing is relatively easy to explain. The total of 16 states at small ω is the maximum number of discrete states the system can have. As ω increases, the coupling between the 16 states strengthens and the smaller basins on the landscape merge into larger basins. As a result, we have less number of states. The optimal transition rate and diversification of transition paths at moderate ω non-adiabatic regime is harder to understand. The unified landscape in extended space picture shed light on this issue. From the 6 dimensional dynamics in extended space we can see we need to consider deterministic dynamics and fluctuations of both protein concentration in x-variables and binding/unbinding processes in c-variables.
At adiabatic limit, the system is close to steady state of H b . Gene states change so frequently that there is not enough time for protein concentration to reach transcription level of the final state. Transitions in x-space are difficult to occur. At small ω non-adiabatic limit, the system is close to steady state of H 0 . System tends to stay in one of the 16 steady states of H 0 and transitions in c-space (change of gene states) are difficult to occur. It is only in moderate ω non-adiabatic regime where H 0 rates and H b rates are comparable, transitions can occur in both x-space and c-space. In this regime transition rate is optimal and transition paths are diversified.