Example 3: Signaling cascade model

The DynGPT package can be applied to solve a given state transition network and infer the associated model parameters. In this section, we will demonstrate the basic use of DynGPT using the signaling cascade model as an illustrative example.

The STN of signaling cascade model can be described as:

1. The state is \(s={{\left( {{s}_{1}},\ldots ,,{{s}_{M}} \right)}^{\text{T}}}\in {{\mathbb{N}}^{\text{M}}}\) with \(M=10,\) where \({{s}_{i}}\) represents the \(i\text{-th}\) signal.

2. Event sets are \(\mathcal{R}=\left\{ \left. {{\mathcal{R}}_{1}},\ldots ,{{\mathcal{R}}_{2M}} \right\} \right.,\) \({{\mathcal{R}}_{1}}\) represents the production the first signal\({{\mathcal{R}}_{i}}\text{ for }i=2,\ldots ,M\) represents the conversion from the current signal to the next signal, \({{\mathcal{R}}_{i}}\text{ for }i=M+1,\ldots ,2M\) represents the signal degradation.

3. State change vectors are: \(\left\{ {{\Delta }^{\left( k \right)}} \right\},k=1,2,\ldots ,20,\) where \({{\Delta }^{\left( k \right)}}=\left( \Delta _{1}^{\left( k \right)},\Delta _{2}^{\left( k \right)},\ldots ,\Delta _{2M}^{\left( k \right)} \right)\) represents the state change vector corresponding to \(k\text{-th}\) event with

   \[\begin{split}\begin{align} & \Delta _{i}^{\left( 1 \right)}={{\delta }_{i,1}}, \\ & \Delta _{i}^{\left( p \right)}=-{{\delta }_{i,p-1}}+{{\delta }_{i,p}}\text{ for }p=2,\ldots ,M, \\ & \Delta _{i}^{\left( q \right)}=-{{\delta }_{i,q-M}}\text{ for }q=M+1,\ldots ,2M, \\ \end{align}\end{split}\]

   where \(i\in \left\{ \text{1},\text{2},\ldots ,2M \right\}\) and \({{\delta }_{i,j}}\) is the Kronecker delta.

4. The propensity functions are: \({{\lambda }_{1}}=\rho ,\text{ }\left( {{\lambda }_{2}},\ldots ,{{\lambda }_{M}} \right)=\left( k{{s}_{1}},\ldots ,k{{s}_{M-1}} \right)\) and \(\left( {{\lambda }_{M+1}},\ldots ,{{\lambda }_{2M}} \right)=\left( \delta {{s}_{1}},\ldots ,\delta {{s}_{M}} \right).\) In this STN, the input of the DynGPT-Solver is the parameters \({{\theta }_{\text{dyn}}}=\left\{ \rho ,k,\delta \right\}\).

1. Generate configuration file

To facilitate data generation and subsequent neural network training, we need create a configuration file for the signaling cascade that encapsulates its key properties. The key propertie of STN primarily consists of the following components:

1. State-related configuration: This includes the names of the states. In signaling cascade model, the state names (state_names) are: [“S1”, “S2”, “S3”, “S4”, “S5”, “S6”, “S7”, “S8”, “S9”, “S10”]

2. Event-related configuration: The occurrence rate functions governing event dynamics, and state transition vectors associated with each event’s occurrence. In signaling cascade, the state transition vectors (state_transition) are: \(\Delta _{i}^{\left( 1 \right)}\), \(\Delta _{i}^{\left( p \right)}\) and \(\Delta _{i}^{\left( q \right)}\). Occurrence rates for each event (event_rates) are: \({{\lambda }_{1}}=r\), \(\left( {{\lambda }_{2}},\ldots,{{\lambda }_{M}} \right)=\left( k{{s}_{1}},\ldots, k{{s}_{M-1}} \right)\), \(\left( {{\lambda }_{M+1}},\ldots ,{{\lambda }_{2M}} \right)=\left( \delta {{s}_{1}}\ldots ,\delta {{s}_{M}} \right)\).

3. Parameter-related configuration: The name of parameter values and the range of parameters. In signaling cascade model,the parameter names (param_names) are: [“ρ”, “k”, “δ”]. Upper bounds of parameters (upper_bound_val) are: [50, 1, 1]. Lower bounds of parameters (lower_bound_val) are: [0.1, 0.1, 0.1].      Additionally, we need to specify several hyperparameters for performing SSA simulations on signaling cascade model, including the model name (model_name), the initial values (initial_val), the maximum simulation time (t_end), the number of parameter sets to simulate for the training (train_data_size) and validation datasets (valid_data_size), and the directory for storing the resulting simulation data (dataset_dir). We define the subset of state indices corresponding to the state names to extract (minimal_state_indices).

Hence, the following configuration files can be defined for signaling cascade model. And we use Python’s json package to store this configuration for subsequent simulation and neural network training.

config_file = {
    "state_names": ["S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9", "S10"],
    "events": [
        {"state_transition": {"S1": 1}, "event_rates": "ρ"},
        {"state_transition": {"S1": -1, "S2": 1}, "event_rates": "k*S1"},
        {"state_transition": {"S2": -1, "S3": 1}, "event_rates": "k*S2"},
        {"state_transition": {"S3": -1, "S4": 1}, "event_rates": "k*S3"},
        {"state_transition": {"S4": -1, "S5": 1}, "event_rates": "k*S4"},
        {"state_transition": {"S5": -1, "S6": 1}, "event_rates": "k*S5"},
        {"state_transition": {"S6": -1, "S7": 1}, "event_rates": "k*S6"},
        {"state_transition": {"S7": -1, "S8": 1}, "event_rates": "k*S7"},
        {"state_transition": {"S8": -1, "S9": 1}, "event_rates": "k*S8"},
        {"state_transition": {"S9": -1, "S10": 1}, "event_rates": "k*S9"},
        {"state_transition": {"S1": -1}, "event_rates": "δ*S1"},
        {"state_transition": {"S2": -1}, "event_rates": "δ*S2"},
        {"state_transition": {"S3": -1}, "event_rates": "δ*S3"},
        {"state_transition": {"S4": -1}, "event_rates": "δ*S4"},
        {"state_transition": {"S5": -1}, "event_rates": "δ*S5"},
        {"state_transition": {"S6": -1}, "event_rates": "δ*S6"},
        {"state_transition": {"S7": -1}, "event_rates": "δ*S7"},
        {"state_transition": {"S8": -1}, "event_rates": "δ*S8"},
        {"state_transition": {"S9": -1}, "event_rates": "δ*S9"},
        {"state_transition": {"S10": -1}, "event_rates": "δ*S10"}
    ],
    "parameters": {
        "param_names": ["ρ", "k", "δ"],
        "upper_bound_val": [50, 1, 1],
        "lower_bound_val": [0.1, 0.1, 0.1]
    },
    "model_name": "signal_cascade",
    "initial_val": [1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    "t_end": 6000,
    "train_data_size": 10000,
    "valid_data_size": 1000,
    "dataset_dir": "/GPUFS/sysu_jjzhang_1/hzw/academicCode/DynGPTTest/DynGPT/data/signal_cascade/",
    "minimal_state_indices": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
    "state_upper_bound": [113, 62, 42, 33, 28, 24, 22, 19, 18, 90]
}
with open("config/{}_config.json".format(config_file["model_name"]), "w") as file:
    json.dump(config_file, file, indent=4)

2. Generate datasets

After generating the configuration file, we construct and simulate a state transition network for the gene expression central dogma model to obtain the dataset required for training the neural network addressing the model. The simulation is implemented in Julia, which provides higher efficiency for stochastic simulation algorithms (SSA) compared to Python. To generate the training data, execute the following command in the terminal:

julia simulation.jl /path/to/config.json

3. Train DynGPT-Solver

# import packages
import sys
import os
import numpy as np
import json
import logging
import warnings
root_dir = "/GPUFS/sysu_jjzhang_1/hzw/academicCode/DynGPTTest/DynGPT/"
sys.path.append(root_dir)
os.chdir(root_dir)
import dyngpt as dg
logging.getLogger().setLevel(logging.ERROR)
warnings.filterwarnings("ignore")

To train a neural network for addressing the signaling cascade model and inferring its parameters, a basic configuration file must first be defined. This file includes four key components: (1) the information of the signaling cascade model, (2) the neural network architecture, (3) the directory for saving outputs and visualizations, and (4) the default parameter settings for the neural network. Subsequent configurations for training, sampling, and inference build upon this foundation by incorporating additional task-specific parameters as needed. The relationships between the configurations of each subsequent stage are presented as follows:

Firstly, the basic configuration file must include detailed information about the dynamic model.

config_basic = {}
config_path_dynmodel = "your_path_config"
model_name ="signal_cascade"
config_path_dynmodel = "config/{}_config.json".format(model_name)
with open(config_path_dynmodel, 'r') as file:
    dynmodel_config = json.load(file)
config_basic = dg.tl.update_dynmodel_config(config_basic,config_path_dynmodel)

Next, the architecture of the neural network is defined by specifying the relevant structural details within the configuration.

config_basic["num_encoder_layers"] = 1  # transformer n_layer
config_basic["n_head"] = 8  # transformer n_head
config_basic["feed_forward_dimension"] = 1024  # transformer feed_forward_dimension

Additionally, we must designate a directory for saving output files and visualizations.

config_basic["dataset_dir"]="your/path/to/saveresult/"
config_basic["figure_dir"]="your/path/to/savefigure/"
config_basic["result_dir"]="result/{}/".format(config_basic["model"])
config_basic["figure_dir"]="figure/{}/".format(config_basic["model"])

Finally, execute the provided code to append default configuration parameters.

config_basic = dg.tl.update_default_config(config_basic)

After defining the configuration file, we proceed to train the neural network designed to solve the signaling cascade. The training process is divided into two stages: the first stage involves pretraining the neural network, while the second stage focuses on fine-tuning it. Additionally, for each stage, we will visualize the corresponding loss function curves to evaluate the training progress.

For model pretraining, the configuration file is extended from the basic configuration. We should specify the number of epochs (epochs_pretrain) and the training dataset path (train_dataset_path) in the configuration file.

config_pretrain = config_basic
config_pretrain.epochs_pretrain = 1000
config_pretrain.train_dataset_path = "your/path/to/train_dataset"
config_pretrain.train_dataset_path = "data/isc/isc_train_stable.json"
result_pre_train = dg.solver.pre_training(config_pretrain)

Visualize the loss values during the pretraining phase of the model.

result_pre_train = np.load("result/signal_cascade/train_loss/signal_cascade_epoch900.npz")
dg.pl.plot_loss(result_pre_train["loss_mean"])

During the fine-tuning phase, the configuration file is also extended from the basic configuration, specifying the number of epochs (epochs_fine_tune), the training dataset path (train_dataset_path), and the pretrained network weights (weight_nn_pretrain_path).

config_fine_tune = config_basic
config_fine_tune.epochs_fine_tune = 1000
config_pretrain.train_dataset_path = "your/path/to/train_dataset"
config_fine_tune.weight_nn_pretrain_path = result_pre_train["weight_nn_pretrain_path"]

result_fine_tune = dg.solver.fine_tune_training(config_fine_tune)

Visualize the loss values during the fine tunning phase of the model.

result_fine_tune = np.load("result/signal_cascade/train_loss/fine_tune_signal_cascade_epoch800.npz")
dg.pl.plot_loss(result_fine_tune["loss_mean"])

After the neural network training converges, we valid the performance of the neural network model in solving signaling cascade from comparative analysis of the samples generated by the neural network and those obtained through simulations in terms of their mean, standard deviation, and overall distribution. When using the trained model for sampling, the configuration file is also extended from the basic configuration, specifying the validation dataset path (valid_dataset_path) and incorporating the fine-tuned network weights (weight_nn_fine_tune_path).

config_sampling = config_basic
valid_dataset_path =  "your/valid/datasetpath/"
valid_dataset_path =  "data/isc/isc_valid_stable.json"
# config_sampling.weight_nn_fine_tune_path = result_fine_tune["weight_nn_fine_tune_path"]
config_sampling.weight_nn_fine_tune_path = "result/signal_cascade/weights/fine_tune_signal_cascade_epoch800.pt"
config_sampling.valid_dataset_path = valid_dataset_path
result_valid_set = dg.tl.compute_sampling_stats(config_sampling)

the indices length is 100

We visualize the performance of the neural network model on the validation set. This includes a comparative analysis of the samples generated by the neural network and those obtained through simulations in terms of their mean, standard deviation, and overall distribution.

dg.pl.plot_model_comparison_stats(result_valid_set,state_index=[0,1,2,3])
dg.pl.plot_model_comparison_stats(result_valid_set,state_index=[4,5,6,7])
dg.pl.plot_model_comparison_stats(result_valid_set,state_index=[8,9])

result_valid_set['observed_datas'],result_valid_set['nn_sample_datas'] = result_valid_set["ssa_sample"],result_valid_set["nn_sample"]
dg.pl.plot_distribution_comparison_nd(result_valid_set)

4. Infer the parameters with DynGPT-Inferrer

In this section, we employ the trained model to infer the corresponding stochastic dynamics from the observed data.

valid_dataset_path =  "data/isc/isc_valid_stable.json"
observed_data,true_params = dg.tl.load_synthetic_data(valid_dataset_path)

For the inference process, the configuration file is also extended from the basic configuration, defining the following hyperparameters: the trained neural network parameters path(weight_nn_fine_tune_path), the loss function (inference_loss_func), the learning rate(inference_lr), the maximum number of iterations(iteration_steps), the initial number of parameters sampled from the prior distribution(random_r0_number), and the threshold for parameter acceptance(param_acceptance_threshold).

config_infer = config_basic
# config_infer.weight_nn_fine_tune_path = result_fine_tune["weight_nn_fine_tune_path"]
config_infer.weight_nn_fine_tune_path = "result/signal_cascade/weights/fine_tune_signal_cascade_epoch800.pt"
config_infer.inference_loss_func = "Likelihood"
config_infer.inference_lr = 0.0001
config_infer.iteration_steps = 80
config_infer.random_r0_number = 2
config_infer.param_acceptance_threshold = 0.5
result_infer = dg.inferrer.inferring_dynamics(observed_data,config_infer,new_model=True,synthetic_flag=True,test_params = 10)

We visualize the distribution obtained from the observational data with the distribution corresponding to the inferred dynamical parameters.

dg.pl.plot_distribution_comparison_nd(result_infer,species_index=[0,1,2,3],data_index = [2])

We visualize the posterior distribution of the inferred dynamical model parameters.

result_infer["param_names"] = np.array(["$ρ$", "$k$", "$δ$"])
dg.pl.plot_param_posterior_dist(result_infer,excluded_indexes=[0],true_values=[],param_indexes=[2,6],width_mm = 140,height_inch_ratio = 0.45)