Dynamical density-functional-theory-based modeling of tissue dynamics: Application to tumor growth

We present a theoretical framework based on an extension of dynamical density-functional theory (DDFT) for describing the structure and dynamics of cells in living tissues and tumors. DDFT is a microscopic statistical mechanical theory for the time evolution of the density distribution of interacting many-particle systems. The theory accounts for cell-pair interactions, different cell types, phenotypes, and cell birth and death processes (including cell division), to provide a biophysically consistent description of processes bridging across the scales, including describing the tissue structure down to the level of the individual cells. Analysis of the model is presented for single-species and two-species cases, the latter aimed at describing competition between tumor and healthy cells. In suitable parameter regimes, model results are consistent with biological observations. Of particular note, divergent tumor growth behavior, mirroring metastatic and benign growth characteristics, are shown to be dependent on the cell-pair-interaction parameters.

Figure 1

Density of the cells (left) and local nutrition concentration (right) over time. We assume that the population growth constant $c_1=1$ and the energy scale in the interaction potential between cells $\beta \varepsilon =1$. The diffusion coefficient ratio $\tilde{D}=1$. The nutrient source is homogeneous with $f\left(\mathbf{r}\right)=1$ and ${\tilde{S}}_n=10$, and the nutrient uptake rate ${\tilde{\lambda }}_n=1$. The area of the domain is $25.6^2$ and $\mathrm{\Delta }x=0.1$.

Figure 2

The average cell density [see Eq. (10)] and the average nutrient density [see Eq. (11)], corresponding to the results in Fig. 1.

Figure 3

The local density of the cells (left three columns, for $\tilde{D}$ = 1, 10, and 100, from left to right) and the nutrient density for $\tilde{D}=1$ (right hand column). The population growth constant $c_1=1$ and the energy scale in the interaction potential between cells $\beta \varepsilon =1$. The nutrient source term has ${\tilde{S}}_n=35$ with $f\left(\mathbf{r}\right)$ given in Eq. (6) and nutrient uptake rate ${\tilde{\lambda }}_n=1$. The area of the system is $25.6^2$, with grid spacing $\mathrm{\Delta }x=0.05$.

Figure 4

The average cell density [see Eq. (10)] and the average nutrient density [see Eq. (11)], corresponding to the results in Fig. 3 when the diffusion coefficient $\tilde{D}$=1, 10, and 100, respectively.

Figure 5

Snapshots of several peak splitting events that occur between the times $t=2.05$ and $t=2.10$. The figures above are in time increments of 0.01 going from top left to bottom right, corresponding to the profiles plotted in the left hand column of Fig. 3, which are for $\tilde{D}=1$.

Figure 6

The linear stability threshold for the two species cells [see Eqs. (31) and (32)] plotted in the total density $\rho \equiv {\rho }_1^b+{\rho }_2^b$ versus concentration $\varphi \equiv {\rho }_1^b/\rho$ plane. The uniform density state is linearly unstable above this line. The top plot shows the curves for $R_{11}=1,R_{22}=1.2,R_{12}=1.1,\beta {\varepsilon }_{11}=\beta {\varepsilon }_{22}=1$ and for varying $\beta {\varepsilon }_{12}$, as given in the key. The curves in the lower plot are for varying $R_{22}=1$, 1.5, 1.7, 1.73, 1.8, and 2. We set the cross-interaction radius $R_{12}=\frac{1}{2}(R_{11}+R_{22})$ and $\beta {\varepsilon }_{11}=\beta {\varepsilon }_{12}=\beta {\varepsilon }_{22}=1$.

Figure 7

Top four panels: plots of $({\rho }_1-{\rho }_2)$, the density profile of the cancer cells minus the density of the healthy cells, at times $t=0.1$, 16, 26, and 34.2. The nutrient uptake rates ${\tilde{\lambda }}_{n1}$=1 and ${\tilde{\lambda }}_{n2}$=1, the population growth constants $c_1=c_2=0.5$, and the threshold nutrient concentration for healthy cells $\alpha =2$. The nutrient source is homogeneous, with $f\left(\mathbf{r}\right)=1$ and ${\tilde{S}}_n=9$. The area of the domain is $25.6\times 25.6$ and $\mathrm{\Delta }x=\mathrm{\Delta }y=0.1$. The cell-cell pair interaction potential parameters are $\beta {\varepsilon }_{11}$=1, $\beta {\varepsilon }_{12}$=1.5, $\beta {\varepsilon }_{11}$=1, $R_{11}=R_{22}=1$, and $R_{12}=0.9$. Bottom: the corresponding average cell density [see Eq. (10)] and the average nutrient density [see Eq. (11)].

Figure 8

Snapshots of $({\rho }_1-{\rho }_2)$, the density profile of the cancer cells minus that of the healthy cells, at the times $t=0.1$, 6.5, 10, and 20 evolving from the initial conditions defined in Eqs. (43) and (44). The system parameters are ${\tilde{\lambda }}_{n1}={\tilde{\lambda }}_{n2}=1,{\tilde{D}}_c={\tilde{D}}_2=1,c_1=c_2=0.5$, and $\alpha =2$. The nutrient source is homogeneous with $f\left(\mathbf{r}\right)=1$ and ${\tilde{S}}_n=9$. The area of the domain is $25.6\times 25.6$ and $\mathrm{\Delta }x=\mathrm{\Delta }y=0.1$. The parameters in the pair interaction potentials between the cells are $\beta {\varepsilon }_{11}=1,\beta {\varepsilon }_{12}=1.5,\beta {\varepsilon }_{11}=1,R_{11}=R_{22}=1$, and $R_{12}=0.9$. In the bottom left panel are plotted the corresponding average cell densities [see Eq. (10)] and the average nutrient density [see Eq. (11)]. In the bottom right panel we plot the trajectory of the time evolution in the $(\rho ,\varphi )$ plane. Note that the points on this trajectory correspond to the integer times $t=0,1,2,\cdots$. We also plot the linear stability threshold for this system. When the trajectory dips below this line, the system temporarily “melts.”

Figure 9

Snapshots of $({\rho }_1-{\rho }_2)$ at the times $t=1.1$, 6.5, 8.5, and 20. All the parameters here are the same as those in Fig. 8, except here the cross interaction pair potential radius is $R_{12}=1$, which is slightly larger.

Figure 10

Snapshots of $({\rho }_1-{\rho }_2)$ at the times $t=1.1$, 5.5, 8.5, and 20. All the parameters here are the same as those in Figs. 8 and 9, except here the cross interaction pair potential radius is even larger, $R_{12}=1.1$.

Figure 11

Snapshots of $({\rho }_1-{\rho }_2)$, for various $\beta {\varepsilon }_{12}=1$ (left), $\beta {\varepsilon }_{12}=1.75$ (middle), and $\beta {\varepsilon }_{12}=2$ (right) and various different times, with time increasing from top to bottom, as indicated above. The other pair potential parameters are $\beta {\varepsilon }_{11}=\beta {\varepsilon }_{22}=1,R_{11}=R_{22}=1$, and $R_{12}=0.9$. The other model parameters are ${\tilde{\lambda }}_{n1}={\tilde{\lambda }}_{n2}=1,c_1=c_2=0.5,\alpha =2$, and ${\tilde{S}}_n=9$ with $f\left(\mathbf{r}\right)=1$. The area of the domain is $25.6\times 25.6$ and $\mathrm{\Delta }x=\mathrm{\Delta }y=0.1$.

Figure 12

On the left are plots of the average cell densities [see Eq. (10)] and the average nutrient density [see Eq. (11)] and on the right plots of the trajectory in the $(\rho ,\varphi )$ plane with the corresponding linear stability threshold line, corresponding to the results in Fig. 11. These are for varying $\beta {\varepsilon }_{12}=1$ (top), $\beta {\varepsilon }_{12}=1.75$ (middle), and $\beta {\varepsilon }_{12}=2$ (bottom).