Simulating gene expression in twin data with a given gene-regulatory network

This tutorial makes a simulated dataset with the structure TwINFER expects: a population of recently divided sister cells, or twins, measured at regular intervals after division. Because we write the network ourselves, the ground truth is known, which makes simulated data the natural way to test the framework.

The simulation runs a population of mother cells until gene expression reaches steady state, divides each mother into a twin pair, and then records both members of every pair over the following 48 hours. The output is a single .csv file, the input table of the inference tutorials (inference_simulation_one_time_point.ipynb and inference_simulation_two_time_points.ipynb).

There are three things to prepare (two input files and one configuration dictionary) and then one function to call.

Step	Section
Choose the network	1. Input file 1: the connectivity matrix
Choose the kinetics	2. Input file 2: the parameter file
Write the settings	3. The configuration dictionary
Simulate	4. Running the simulation
Check the result	5. The output and 6. Confirming steady state

Note. Which time points are compared (one time point, or t₁ and t₂) and how the twins are split into groups are chosen at inference, with time_points in infer_with_twinfer. The simulation records every cell at every time point, so the timing decisions here are to simulate long enough after division to cover the times we intend to use later and the resolution of measurement.

Setup

Install the package once from the repository root (pip install -e package/, see the README) and start Jupyter from the repository. Then:

Note. The first simulation of a session takes an extra minute or two while numba compiles the simulation kernels. Later runs are faster.

import json
import os

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from twinfer.utils.paths import get_repo_root

path_to_code_repo = get_repo_root()
path_to_data = path_to_code_repo / "simulation_example_input_data"
sorted(p.name for p in path_to_data.iterdir())

['connectivity_matrix_2_gene_no_interaction.txt',
 'connectivity_matrix_3_gene_linear_cascade.txt',
 'connectivity_matrix_4_gene_linear_cascade.txt',
 'connectivity_matrix_5_gene_linear_cascade.txt',
 'connectivity_matrix_A_B.txt',
 'connectivity_matrix_A_and_B.txt',
 'connectivity_matrix_A_and_B_both_repress.txt',
 'connectivity_matrix_A_rep_B.txt',
 'connectivity_matrix_A_represses_B_B_activates_A.txt',
 'connectivity_matrix_A_to_B.txt',
 'connectivity_matrix_fan_in.txt',
 'connectivity_matrix_fan_out.txt',
 'connectivity_matrix_feed_forward.txt',
 'connectivity_matrix_figure_1_network.txt',
 'connectivity_matrix_mutual_regulation.txt',
 'median_parameter.csv']

1. Input file 1: the connectivity matrix

Action. Choose one of the matrices in simulation_example_input_data/, or write your own. The configuration below points at connectivity_matrix_A_to_B.txt, so this step can be skipped on a first run.

A plain text file with the network topology: a square grid of integers, one row and one column per gene, separated by commas. Take row i and column j; the entry is the effect of gene i on gene j.

Value	Meaning
`0`	no interaction
`1`	gene i activates gene j
`-1`	gene i represses gene j

The diagonal is always 0; self-regulation is not represented in this model.

Important. Rows are the regulator and columns are the target. Transposing the matrix reverses every interaction in the network, so check this before starting a long simulation: a transposed file gives a perfectly valid simulation of the wrong network.

A two-gene example

The complete contents of connectivity_matrix_A_to_B.txt:

print((path_to_data / "connectivity_matrix_A_to_B.txt").read_text())

Reading this row by row:

Row 1 is 0,1. The first entry is the effect of gene 1 on itself, which is always 0. The second entry is the effect of gene 1 on gene 2, and it is 1, so gene 1 activates gene 2.
Row 2 is 0,0. Gene 2 has no effect on gene 1, and none on itself. There is no reverse interaction.

The network is therefore a single activating edge, gene 1 → gene 2. To make it repressive instead, change the 1 to -1; this is the file connectivity_matrix_A_rep_B.txt:

print((path_to_data / "connectivity_matrix_A_rep_B.txt").read_text())

A three-gene example

connectivity_matrix_feed_forward.txt:

print((path_to_data / "connectivity_matrix_feed_forward.txt").read_text())

Row 1: gene 1 activates gene 2 and gene 3.
Row 2: gene 2 activates gene 3.
Row 3: gene 3 regulates nothing.

This is the feed-forward loop, one of the motifs that correlation-based methods commonly misclassify.

Other matrices provided

simulation_example_input_data/ also contains matrices for the no-interaction control, bidirectional and mixed-sign two-gene networks, three- to five-gene cascades, fan-in and fan-out motifs, mutual regulation, and the network of Figure 1 (listed in the Setup cell above). The names are descriptive, and each file can be read with the rules above.

Writing a new matrix

Write a grid of 0, 1 and -1 separated by commas, one matrix row per line, and save it with a .txt extension.

The easiest route is to start from an existing file rather than a blank one: copy whichever of the provided matrices is closest in size to the network you want, open it in a text editor, and edit the entries. This keeps the format correct and only the topology changes.

2. Input file 2: the parameter file

Action. Nothing, on a first run. median_parameter.csv is supplied with the code and is used as it is; the only decision is which of its rows to give each gene, set through rows_to_use in the configuration below.

The connectivity matrix says which interactions exist. The parameter file says how strong they are, and how each gene behaves without regulation.

Each gene needs six parameters, covering transcriptional bursting, transcription, translation, and the stability of mRNA and protein. Each regulatory interaction needs two more.

The file holds these as a table. Every row is one complete parameter set, and the rows are numbered 0, 1, 2, and so on. We do not edit the file to choose a parameter set; rows_to_use in the configuration names the row each gene takes.

For most purposes median_parameter.csv can be used unchanged. Its values are the medians and defaults compiled from published measurements in mammalian cells.

parameters = pd.read_csv(path_to_data / "median_parameter.csv", index_col=0)
parameters

	k_on	k_off	mrna_half_life	protein_half_life	k_prod_mRNA	k_prod_protein	n_gene_1_to_gene_2	k_add_gene_1_to_gene_2	n_gene_2_to_gene_1	k_add_gene_2_to_gene_1	n_gene_1_to_gene_3	k_add_gene_1_to_gene_3	n_gene_2_to_gene_3	k_add_gene_2_to_gene_3	n_gene_3_to_gene_2	k_add_gene_3_to_gene_2
0	0.66	8.6	4	45	2	560	2	6.0	2	6.0	2	6.0	2	6.0	2	6.0
1	0.66	8.6	4	45	2	560	2	6.0	2	6.0	2	6.0	2	6.0	2	6.0
2	0.12	8.6	4	45	2	560	2	6.0	2	6.0	2	6.0	2	6.0	2	6.0
3	1.66	8.6	4	45	2	560	2	6.0	2	6.0	2	6.0	2	6.0	2	6.0
4	0.66	8.6	4	45	2	560	2	0.8	2	0.8	2	0.8	2	0.8	2	0.8
5	0.66	8.6	4	45	2	560	2	0.8	2	6.0	2	0.8	2	0.8	2	0.8
6	0.66	8.6	4	45	2	560	2	2.0	2	2.0	2	2.0	2	2.0	2	2.0

Column	Description
`k_on`	Rate of burst initiation, per hour. Median 0.66. The gene switches stochastically between an inactive and a short-lived active state, and is transcribed only in the active state.
`k_off`	Rate of burst inactivation, per hour. Median 8.6.
`mrna_half_life`	mRNA half-life in hours. Median 4. The degradation rate is derived from it assuming first-order decay.
`protein_half_life`	Protein half-life in hours. Median 45.
`k_prod_mRNA`	Transcription rate while the gene is active, in mRNA per hour. Median 2.
`k_prod_protein`	Translation rate, in protein per mRNA per hour. Median 560.
`n_gene_i_to_gene_j`	Hill coefficient of the interaction, setting how sharply the response saturates. Default 2.
`k_add_gene_i_to_gene_j`	Maximum change in burst frequency produced by the interaction, in h⁻¹. The defaults are 6 for activation and 0.8 for repression. The sign of the interaction comes from the connectivity matrix, not from this value, which is always positive.

Which row to use

Row 0 and row 4 of median_parameter.csv are identical except for k_add:

row	`k_on`	`k_add`	intended use
0	0.66	6	activating networks
4	0.66	0.8	repressing networks

So a two-gene activation network uses 'rows_to_use': [[0, 0]], and the corresponding repression network uses [[4, 4]]. Rows 2 and 3 differ from row 0 only in k_on (0.12 and 1.66), which is how populations with two distinct transcriptional states are built.

Two further points affect how the file is read:

The Hill constant is not taken from this file. It is set automatically to the steady-state transcription factor level of the regulator.
If a k_add column is missing for some edge, the sign-based default (6 or 0.8) is applied automatically. Setting use_csv_k_add to False in base_config applies the defaults to every edge.

Important. Row indices start at 0, and each row must have a distinct index, since the index appears in the output filename.

3. The configuration dictionary

Action. Adjust four things in the cell below: the file in path_to_connectivity_matrix, rows_to_use so that it has one value per gene, type so the output files are recognisable, and path_to_output, the folder the simulation is written to. Everything else can stay at the values shown for a first run.

All settings live in one dictionary, base_config.

# Folder the simulation file and its log are written to (created if needed).
path_to_output = path_to_code_repo / "simulation_example_output_data"

base_config = {
    # ---- population size and timing ----
    'n_cells': 6000,
    'simulation_time_before_division': 1500,
    'twin_simulation_time_after_division': 48,
    'twin_measurement_resolution': 1,

    # ---- network and parameters ----
    'path_to_connectivity_matrix': str(path_to_data / "connectivity_matrix_A_to_B.txt"),
    'param_csv': str(path_to_data / "median_parameter.csv"),
    'rows_to_use': [[0, 0]],
    'combinatorial_interaction_type': "additive",

    # ---- naming and file output ----
    'type': "A_to_B",
    'output_folder': str(path_to_output),
    'log_file': str(path_to_output / "log.jsonl"),
    'save_before_division': False,

    # ---- computational resources ----
    'number_of_parallel_parameters': 1,
    'number_of_cores_per_parameter': 4,
}

Settings reference

Setting	Description
`n_cells`	Number of mother cells, and therefore the number of twin pairs. The output contains twice this many cells. The default of 6000 pairs parallels typical scRNA-seq capabilities and is the number of twin pairs used in the paper’s simulations.
`simulation_time_before_division`	Simulated hours for which the mother cells evolve before dividing. This brings the population to steady state and is not recorded. With the median parameters a thousand hours or more is required; 1500 is a safe default. Increase it if the steady-state check (section 6) fails.
`twin_simulation_time_after_division`	Simulated hours over which the twins are followed after division, and the period that appears in the output. It must cover the latest time point you intend to use at inference. A value of 48 is comfortable.
`twin_measurement_resolution`	Interval between recorded measurements, in hours. A value of `1` gives 49 measurements per cell over 48 hours, from hour 0 through hour 48. Finer resolution gives proportionally larger files.
`path_to_connectivity_matrix`	Full path to the connectivity matrix.
`param_csv`	Full path to the parameter file.
`rows_to_use`	Assignment of parameter rows to genes. `[[0, 0]]` specifies one simulation in which both genes take the parameters of row 0. The inner list must contain one entry per gene, so a three-gene network requires `[[0, 0, 0]]`. Additional inner lists specify additional simulations: `[[0, 0], [4, 4]]` runs two.
`combinatorial_interaction_type`	How several regulators of the same target are combined. `"additive"` sums the pairwise contributions, `"AND"` requires all regulators to be bound, `"OR"` requires any one. Must be spelled exactly `"additive"`, `"AND"` or `"OR"` (anything else raises an error). It has no effect when no gene has more than one regulator.
`type`	Short identifier for the network, used in the output filenames. Avoid spaces.
`output_folder`	Directory in which the output `.csv` is written.
`log_file`	Path to a running log in JSON Lines format. One entry is appended per simulation, recording the parameters used.
`save_before_division`	Whether to also write the trajectories of the mother cells before division (`simulation_before_division_df_*.csv`, every cell at every hour of `simulation_time_before_division`). This file is large and is not used by the inference; it is switched off here. The default, if the key is left out, is `True`.
`number_of_parallel_parameters`	Number of parameter sets run at the same time. Relevant only when `rows_to_use` contains more than one entry.
`number_of_cores_per_parameter`	CPU cores given to each simulation. This value multiplied by `number_of_parallel_parameters` should equal the number of cores available.

Warning. Key names in base_config are not validated. A misspelled key is ignored without an error, and the simulation silently uses the default value instead. For example, a misspelled combinatorial_interaction_type key makes a network meant to use "AND" or "OR" logic run as "additive", and nothing in the output shows that the setting was dropped. Check the spelling against the table above.

Optional settings

These keys are not required; leaving them out gives the standard behaviour.

Setting	Default	Description
`divide_binomial`	`False`	If `True`, transcripts and proteins are partitioned binomially between the two daughter cells at division rather than perfectly, so the twins differ slightly at time zero.
`p_major`	`0.5`	Probability that a given transcript or protein goes to the first daughter cell. `0.5` is symmetric partitioning; larger values give asymmetric partitioning, which lies outside the assumptions of the framework.
`use_csv_k_add`	`True`	Whether `k_add` values are read from the parameter file. `False` applies the sign-based defaults to every edge.
`k_add_matrix`, `n_matrix`	`None`	n × n arrays of per-edge interaction strengths and Hill coefficients, used for edges absent from the parameter file.
`use_given_K`, `K_to_use`	`False`, `None`	Supply the Hill constants directly as an n × n matrix instead of deriving them from the steady state.
`scale_K`	`None`	An n × n array of multipliers applied to the derived Hill constants.

4. Running the simulation

Action. Run the two cells below in order. The second one is the simulation itself and is the slow step. The simulation takes about 15 minutes for 6000 parent cells with 1500h to achieve steady state using 4 cores.

Set the number of threads and load the simulation function:

from numba import get_num_threads, set_num_threads

set_num_threads(base_config['number_of_cores_per_parameter'])
print("Threads Numba will use:", get_num_threads())

from twinfer.simulation.gillespie_simulations import process_param_set

Threads Numba will use: 4

Note. The thread count can be set only once per session. If number_of_cores_per_parameter is changed afterwards, restart the kernel before running again.

Run it:

os.makedirs(base_config['output_folder'], exist_ok=True)

rows_to_use = base_config['rows_to_use']
labels = ["rows_" + "_".join(map(str, row)) for row in rows_to_use]

path_to_simulation_file = process_param_set(rows_to_use[0], labels[0], base_config)
print(f"Saved the simulation file as {path_to_simulation_file}")

Log pi on is set to False
Combinatorial interaction type: additive

k_add resolution summary
Target              Scope                   Source         Value
----------------------------------------------------------------
gene_2              per-gene                CSV avg       6.0000

Removed 5 unused k_add key(s): ['{k_add_gene_1_to_gene_2}', '{k_add_gene_1_to_gene_3}', '{k_add_gene_2_to_gene_1}', '{k_add_gene_2_to_gene_3}', '{k_add_gene_3_to_gene_2}']

Starting base simulation

process_param_set takes one parameter-row assignment ([0, 0]), a short label used in the filename (built automatically as "rows_0_0"), and the configuration dictionary. It returns the path of the file it wrote, which is the input of the inference tutorials.

While it runs, it prints a steady-state check of the mother cells at the end of the pre-division period. If the population does not pass, it prints a warning and adds an entry to error_log.jsonl in the folder of the log file; see section 6.

To run every parameter set of rows_to_use at the same time instead of only the first:

from joblib import Parallel, delayed

paths = Parallel(n_jobs=base_config['number_of_parallel_parameters'])(
    delayed(process_param_set)(row, label, base_config)
    for row, label in zip(rows_to_use, labels)
)

5. The output

Action. Nothing needs to be configured here. Keep the path printed above: it is the input of the inference tutorials.

Filenames

One .csv is written per parameter set, with an automatically generated name, for example:

df_rows_0_0_24102026_140638_ncells_6000_A_to_B_04cf9c46.csv
   │       │              │           │       └── identifier of this run
   │       │              │           └────────── the 'type' identifier
   │       │              └────────────────────── number of twin pairs (n_cells)
   │       └───────────────────────────────────── date and time (DDMMYYYY_HHMMSS)
   └───────────────────────────────────────────── parameter rows used

The identifier is random for every run, so runs with different settings cannot overwrite one another. The same identifier appears in log.jsonl, which allows the parameters behind any output file to be recovered.

What the output file contains

df = pd.read_csv(path_to_simulation_file)
print(df.shape)
df.head()

Each row records one cell at one measurement time.

Column	Description
`clone_id`	Identifies the twin pair. Two cells sharing a `clone_id` descend from the same mother cell. Runs from 0 to `n_cells` − 1.
`replicate`	Which member of the pair the row refers to, `1` or `2`. Together with `clone_id` this defines the twin relationship the framework depends on.
`cell_id`	Unique identifier per cell, running from 0 to 2 × `n_cells` − 1.
`time_step`	Hours since division; `0` is the moment of division.
`gene_i_mRNA`	Transcript count of gene i, one column per gene. Inference is performed on these counts.
`gene_i_protein`	Transcription factor count of gene i. It drives the regulation during the simulation and is used by the steady-state check.
`gene_i_A`, `gene_i_I`	Promoter state of gene i: `gene_i_A` is 1 while the gene is active (transcribed) and `gene_i_I` is 1 while it is inactive.

The two members of the first twin pair, at division:

df[(df['clone_id'] == 0) & (df['time_step'] == 0)]

Both rows have clone_id 0, with replicate 1 and 2 (cell_id 0 and cell_id n_cells). At time_step 0 their counts are identical, because partitioning is perfect by default; their trajectories diverge afterwards, and it is this gradual divergence that TwINFER uses.

With the default settings the file contains 588,000 rows (12,000 cells at 49 measurement times).

What the log file contains

log.jsonl receives one line per simulation, with the run identifier, parameter rows, number of twin pairs, network identifier, timestamp, every rate constant used (including the derived degradation rates and Hill constants) and the predicted steady-state transcription factor (protein) levels. The entry of the run above is the last line:

with open(base_config['log_file']) as f:
    log_entry = json.loads(f.readlines()[-1])

print({key: log_entry[key] for key in ["id", "rows", "n_cells", "type", "timestamp", "steady_state"]})
log_entry["param_dict"]

6. Confirming steady state

Action. Run the cell below and look at the plot before moving on to inference.

The framework assumes the population had reached steady state at division. If it had not, mean expression drifts over the measurement window and the inference will be unreliable. Plotting mean expression against time is a direct check: the curves should be flat. The dashed lines are the predicted steady-state levels from the log file, an analytical approximation given for comparison; the simulated means can sit a few percent away from them.

protein_cols = [c for c in df.columns if c.endswith("_protein")]
mean_protein = df.groupby("time_step")[protein_cols].mean()

fig, ax = plt.subplots(figsize=(5, 3.5))
for col, predicted in zip(protein_cols, log_entry["steady_state"]):
    line, = ax.plot(mean_protein.index, mean_protein[col], label=col)
    ax.axhline(predicted, color=line.get_color(), linestyle="--", linewidth=1)
ax.set_xlabel("Time after division [h]")
ax.set_ylabel("Mean protein per cell")
ax.legend(frameon=False)
plt.show()

pd.DataFrame({"mean at 0 h": mean_protein.iloc[0], "mean at last time step": mean_protein.iloc[-1],
              "predicted steady state": log_entry["steady_state"]})

If the curves trend upward or downward, increase simulation_time_before_division (for example to 3000) and simulate again.

Troubleshooting

“The number of parameter rows entered is less than the number of genes.” rows_to_use contains fewer entries than the network has genes. A three-gene network requires [[0, 0, 0]] rather than [[0, 0]].

The simulation is slow. Run time scales with n_cells and with simulation_time_before_division. To confirm a configuration is valid before a full run, reduce n_cells to 500 for a trial, then restore it.

The output file cannot be found. process_param_set both returns and prints the full path of the file it wrote.

Repressive interactions appear to have no effect. Check that the entry in the connectivity matrix is -1 rather than 1, that it is in the intended position (regulator on the row, target on the column), and that the parameter row gives a k_add meant for repression (row 4). The effective burst frequency is floored at zero, so an unusually large repression strength silences the target rather than producing a stronger signal.

Checklist

Step	Input	Output
Install the package	`pip install -e package/`	`import twinfer` works.
Choose a connectivity matrix	A `.txt` grid of `0`, `1`, `-1`; rows are regulators, columns are targets	The network to simulate.
Choose a parameter file	`median_parameter.csv`, unchanged in almost all cases	The kinetics of each gene and interaction.
Define `base_config`	Paths, `rows_to_use` with one entry per gene, `type`	A dictionary. Nothing is printed.
Set the thread count and import `process_param_set`	`number_of_cores_per_parameter`	`Threads Numba will use: 2`. Restart the kernel first if the core count has changed.
Run `process_param_set`	One entry of `rows_to_use`, a label, and `base_config`	One `.csv` in `output_folder`, one line added to `log.jsonl`, and the file path printed and returned.
Plot mean expression against time	The simulation `.csv`	One line per gene. Flat curves confirm steady state; a trend means `simulation_time_before_division` should be increased.

Next

Score the simulated gene pair with TwINFER: inference_simulation_one_time_point.ipynb (one time point) and inference_simulation_two_time_points.ipynb (two time points, direction of regulation). Both read the example simulation in simulation_example_output_data/; to use the file made here, read path_to_simulation_file instead of the example file in their setup cell.