Examples

Signal Cascade Markov Model

In this example we want to solve a Markovian Masterequation corresponding to a genetic signal cascade. We will use an implicit Euler
method in the time direction using the ALS algorithm to solve the individual steps. The construction of the operator will be done
according to the SLIM decomposition derived in [P. Gelß et al., 2017] (cf. Example 4.1 therein).

Transition Matrices

Our solution tensor $X[i_1, i_2, \dots]$ represent the likelyhood, that there are $i_1$ copies of protein one, $i_2$ copies of protein two
and so on. As the likelyhood of very large $i_j$ becomes small very fast we can restrict ourselves to a finite tensor with $i_j \in \{0,1,\dots,n_j\}$
represented in our sourcecode as

C++Python

constsize_tMAX_NUM_PER_SITE=32;

MAX_NUM_PER_SITE=32

For the different events we can now describe matrices that have the corresponding action. We will denote as $M$ the creation of a
new protein (remove current number of proteins = diagonal equals -1; add current number + 1 proteins = offdiagonal equals 1)

The probability of construction of a protein $x_i$ is actually given in terms of the number of proteins $x_{i-1}$ as $\frac{x_{i-1}}{5+x_{i-1}}$,
so we will need another matrix $L$ that gives these probabilities, such that we can later construct the corresponding two-site TT Operator as $L\otimes M$.

The corresponding destruction could be expressed similarly, but as the probability of destruction in our example only depends on the
number of proteins $x_i$ themselves, we will use a matrix denoted as $S$ instead, that already includes these propabilities.
Relative to the number of proteins $x_i$ the destruction probability in this example can be given as $0.07 x_i$, we thus have:

defcreate_S():S=xe.Tensor([MAX_NUM_PER_SITE,MAX_NUM_PER_SITE])# set diagonalforiinxrange(MAX_NUM_PER_SITE):S[[i,i]]=-i# set offdiagonalforiinxrange(MAX_NUM_PER_SITE-1):S[[i,i+1]]=i+1return0.07*S

System Operator

The operator corresponding to the full system of $N$ proteins can now be expressed with the use of these matrices. As can be
seen in above mentioned paper, it is given by

where includes the construction of the first protein (that does not depend on any other proteins) as $S^* = 0.7\cdot M + S$.

The construction of this operator in xerus is straight-forward: we first construct the individual matrices, use them to
construct the components as given in the previous formula and then simply set them via .set_component as the components of
our operator.

Implicit Euler

To solve the Masterequation in the time domain we will use a simple implicit Euler method. In every step we have to solve
$ (I-\tau A) x_{i+1} = x_i $ for $x_{i+1}$. To do so, we will use the xerus builtin ALS method. From previous experiments
we know, that the _SPD (for “symmetric positive semi-definite”) variant of the ALS works fine in this setting even though
the operator is not symmetric. To define the parameters only once, we create our own ALS variation.

We will keep an eye on the residual after each step for the purpose of this example to ensure, that these claims actually
hold true, and will store the result of every step to be able to plot the mean concentrations over time in the end.

To ensure, that the entries of the solution tensor actually represent probabilities we will also normalize the tensor at every
step to ensure that its one-norm is equal to 1. This norm is usually hard to calculate, but under the assumption that all entries
are positive we can express it as a simple contraction with a ones-tensor.

All we have to do now, is to provide a starting point (no proteins with probability 1) and call the appropriate functions.
As our implicit Euler method as it stands has no means of adapting the rank of the solution we will increaseit from the start
by adding a small pertubation to the starting configuration. To do this we have to convert the (sparse) dirac TTTensor to
dense representation because sparse TT cannot (yet) be added to other TTTensors (see the TTTensor documentation).

Output

At this point we have calculated the solution tensors and can start ot calculate quantities of interest from them. For the
purpose of this example we will simple calculate the mean concentration of every protein at every timestep. To calculate the
mean we simply have to weight the mode corresponding to the protein in question with the number of proteins it represents
(the vector $(0, 1, 2, \dots)$) and sum over all other protein configurations (ie. contract a ones-vector to those modes).
We will use the most general TensorNetwork class to write these contractions. This way xerus can be very lazy and only
perform any actions (and decide upon a contraction order) when we query it for the final value.

C++Python

doubleget_mean_concentration(constTTTensor&_res,constsize_t_i){constIndexk,l;TensorNetworkresult(_res);constTensorweights({MAX_NUM_PER_SITE},[](constsize_t_k){returndouble(_k);});constTensorones=Tensor::ones({MAX_NUM_PER_SITE});for(size_tj=0;j<_res.degree();++j){if(j==_i){result(l&0)=result(k,l&1)*weights(k);}else{result(l&0)=result(k,l&1)*ones(k);}}// at this point the degree of 'result' is 0, so there is only one entry
returnresult[{}];}

defget_mean_concentration(x,i):k,l=xe.indices(2)result=xe.TensorNetwork(x)weights=xe.Tensor.from_function([MAX_NUM_PER_SITE],lambdaidx:idx[0])ones=xe.Tensor.ones([MAX_NUM_PER_SITE])forjinxrange(x.degree()):ifj==i:result(l&0)<<result(k,l&1)*weights(k)else:result(l&0)<<result(k,l&1)*ones(k)# at this point the degree of 'result' is 0, so there is only one entryreturnresult[[]]

Observing the evolution of concentrations over time is now a simple matter of iterating over all solution steps and proteins.