Simulation Process
How to build up a simulation to simulate the profile of a tumor?
The difficulty in simulating an avascular tumor cell is: how to simulate a system with the boundary condition fixed at infinity $u(\pm \infty) = C_0$?It is impossible to create a parameter space with infinite size. We need to transform the original parameter space into a finite space to overcome this problem. In our work, we choose a hyperbolic transformation. That is assuming the original space is denoted byy $x$, and the transformed space $\eta$, say eta-space, would be $$ \eta = \tanh\left(\frac{x}{x_0}\right). $$ By doing so, the boundary condition will turn out to be $u(x = \pm \infty) = u(\eta = \pm 1) = C_0$. However, since hyperbolic tangent is not a linear function, it is necessary to be cautious when we conduct the simulation. To be more specific, when the system moves from $\eta$ to $\eta+0.1$, the trajectory in real space moves from $x_0\text{arctanh}(\eta)$ to $x_0\text{arctanh}(\eta+0.1)$. For $\eta = 0.1$ and $0.6$, the system moves $0.1x_0$ and $0.17x_0$, respectively. Apparently, as the system moves far from the origin when they cross a unit step in eta space, they cross a huge step in real space. To trap the system near the origin such that the crossing step is linear between real space and eta space, we set $x_0 = 400$ in our simulation. See Fig. 1.
Now, we know how to simulate in eta-space.However, how do we find a spatial derivative precisely? Remind that if the system is periodic, the best way to simulate is using Fourier transformation and calculating the derivative in Fourier space. What is the case in the current simulation? Fortunately, the answer appears if we transform the laplacian operator into the eta space. That is, $$ \begin{align*} \nabla^2 &= \frac{\partial }{\partial x}\frac{\partial }{\partial x} = \frac{\partial }{\partial x}\left(\frac{\partial \eta}{\partial x}\frac{\partial }{\partial \eta}\right) \\[1em] &= \frac{\partial }{\partial x}\left(\frac{\text{sech}^2(x/x_0)}{x_0}\frac{\partial }{\partial \eta}\right) \\[1em] &= \frac{\partial }{\partial x}\left(\frac{1-\eta^2}{x_0}\frac{\partial }{\partial \eta}\right) \\[1em] &= \left(\frac{(1-\eta^2)}{x_0^2}\frac{\partial }{\partial \eta}\right)\left((1-\eta^2)\frac{\partial }{\partial \eta}\right) \\[1em] &= \frac{(1-\eta^2)}{x_0^2}\hat{\cal{L}}, \end{align*} $$ where $\hat{\cal{L}}$ is the Legendre operator. Apparently, the best way to find spatial derivative is in Legendre space in this case. Therefore, the most important thing in this simulation is finding the Legendre coefficients of all functions. To get the solution as precise as possible, we utilize the orthonormal property of Legendre polynomial $P_m(x)$, $$ \int_{-1}^1 P_m(x)P_n(x)dx = \frac{2}{2n+1}\delta_{mn} $$ and employ the Gaussian-Legendre quadrature to find the integration. See Wikipedia for more details. In our simulations, we split the space into $4096$ data points, which means we need to process a $4096\times 4096$ matrix multiplication in each simulation step. The size of the matrix suggests we use GPU simulation to accelerate. In the following, I will present some important parts of my program and briefly discuss the codes.
Presentation of my code
First of all, to simulate a function precisely, it would be better to transform between different spaces as less as possible. Therefore, since we need the information of the function at the site of the root of the Legendre polynomial when we employ Gaussian-Legendre quadrature, we simulate all things on the site of the root of the Legendre polynomial, which means the initial conditions of tumor cell density and nutrient density is set on the site of the root of Legendre polynomial. Fortunately, Numpy provides a package for us to find the root of the Legendre polynomial. The code of finding the root and the weighting function $w$ is (Notice that finding $w$ and root is time-consuming if you do not change the split number of the space, I suggest you save the arrays of $w$ and roots):
def findlagendreweight(split):
if os.path.isfile('root_{}.npy'.format(split)):
Broots = np.load('roots_{}.npy'.format(split))
else:
A = np.zeros(split+1)
A[-1] = 1
B = npl.Legendre(A)
Broots = B.roots()
np.save('roots_{}.npy'.format(split), Broots)
# now and after are used to find the derivative of legendre polynomial
now = np.zeros(split+1)
now[-1] = 1
after = np.zeros(split+1)
after[-2] = 1
pder = split*(Broots*npl.legval(Broots, now)-npl.legval(Broots, after))/(Broots**2 -1)
warray = 2/((1-Broots**2)*(pder**2))
np.save('warray_{}.npy'.format(split), warray)
return warray
Next, since we need to multiply the Legendre polynomial to find the coefficient, we calculate all polynomial values in a matrix and save it. It is worth noting that this process is the most time-consuming in this simulation. I strongly suggest you save the array first.
legendrearray = np.zeros(split+1)
legendrearray = []
for i in range(0, split+1):
nowcoeff = np.zeros(i+1)
nowcoeff[-1] = 1
legendrearray.append(npl.legval(xaxis, nowcoeff))
legendrearray = np.array(legendrearray)
np.save('./legendrearray_{}.npy'.format(split), legendrearray)
Then, you can calculate the derivative and use the fourth-order Runge-Kutta method to simulate. The result should be something like Fig. 2
However, as you see, the system always moves out. Since the simulation should not be taken far from the origin, we trap the peak of the tumor cell density profile and calculate its velocity. By transforming the system into a moving frame with the opposite direction of the velocity, the tumor profile will stay stationary (One can imagine it as a treadmill ). There are two points we should be cautious about:
- The moving frame should be used after two surfaces are generated to avoid interactions between them.
- The moving frame has opposite directions in positive ηη and −η−η since the tumor moves along different directions in these two regions.
To track the position of the peak of the tumor profile, we calculate the position for the zero point of the spatial derivative of the tumor profile. Since we know the Legendre coefficients, we can calculate the value of the function anywhere. The process uses the Newton method to find the zero point of $\partial n/\partial x= 0$. See Wikipedia for more details about the Newton method.
def findmax( function, coeff, eigvalue1, eigvalue2, startx = 0):
x = startx
alson = coeff*eigvalue1
nminuscoeff = np.zeros(len(eigvalue1))
nminuscoeff[0:len(nminuscoeff)-1] = coeff[1:len(coeff)]*eigvalue1[1:len(eigvalue1)]
while True:
single = npl.legval(x, alson)*x/(x**2 -1) + npl.legval(x, nminuscoeff)*(-1)/(x**2-1)
if abs(single) < 10**(-12):
break
double=single*2*x/(1-x**2) + npl.legval(x, coeff*eigvalue2)/(1-x**2)
dx=single/double
x -=dx
return x
Finally, if you simulate the system until the profile is unchanged, you get a steady solution of the system. Or, you can utilize the nutrient evolution equation to help double-check the profile. That is, for $\partial u/ \partial t=0$, $$ n = \frac{D_u \nabla^2 u }{bA(u)}. $$ If the simulation profile of $u$ and $n$ satisfy the equality, we can argue that the system has reached its steady state.
Then, if you are interested in finding the instability of the perturbation, you can apply a perturbation with wave number $k$ in $y$ direction to the steady state. For simplicity, I used the $10^{-6}$ times the steady state as the initial state of the perturbation $n_1(t=0)$. By evolving the perturbation until $n_1(t=t)/n_1(t=t+dt)$ reaches the same constant for non-zero $n$ at each site, we can argue that profile is shape preserving. The shape-preserving perturbation is shown in Fig. 3.
Finally, by tracking the peak of the shape-preserving profile, you can know if it grows or decays and judge the instability of the steady state. By conducting a lot of wave numbers, you can know the stability of a steady tumor at that boundary condition. By changing the boundary condition, you can plot all the phase diagrams yourself!