How to run simulations

Monte Carlo simulations are set up, controlled and run through the MC class. There is a juypter notebook in notebooks/tutorial.ipynb that serves as a tutorial on how to get started with the MC class.

Attributes and their meaning

The most important variables to set and their meaning are:

domain: Domain to run the simulation in. Can be 1D, 2D or 3D. Should be a uint array. Labels of this array are used to look up the diffusivities and permeabilities (see below). Label 0 is reserved for the background, which means there cannot be any particles there.
instance_segmentation: Instance segmentation of the domain. Should be a uint array. Adds option to have barriers within the same compartment (e.g. cell membranes).
diffusivities: Diffusivities of the compartments in the domain. 1D array of length n_compartments. Should be in standard units. Labels of the domain are mapped to these diffusivity values (e.g. label 1 --> diffusivities[1]).
permeability: Permeabilities between different compartment boundaries. Needs to be a square matrix of size n_compartments x n_compartments. Should be in standard units. Rows are the source compartments and columns are the target compartments. E.g. permeability[1, 2] is the permeability from compartment 1 to compartment 2. Diagonal elements will be checked if there is a particle crossing from one instance to another within the same compartment (e.g. cell membrane). Permeability values from and to compartment 0 should be set to 0.
permeation_prob: Permeation probabilities of the compartments. Will be calculated from the permeability matrix. Should be a square matrix of size n_compartments x n_compartments. Rows are the source compartments and columns are the target compartments. E.g. permeation_prob[1, 2] is the permeation probability from compartment 1 to compartment 2.
T2: T2 relaxation times of the compartments. 1D array of length n_compartments. Should be in standard units. Labels of the domain are mapped to these T2 values (e.g. label 1 --> T2[1]).
n_walkers: Number of walkers to simulate. Should be an integer.
dt: Time taken for a single particle step. Should be in standard units. This will get converted to a steplength based on the diffusivity and voxel size. dt needs to be small enough to not get a bigger steplength than the voxel size.
voxel_size: Size of the voxels in the domain. Should be in standard units. Currently expects isotropic voxels.
do_reconstruction: If True, the signal will be reconstructed from the walkers (see section Reconstruction).
reconstruct_downsample: Factor to downsample the signal for reconstruction. E.g. if value is 10 and original domain size was 100x90x80 the downsampled domain will be 10x9x8.

How to run a simulation

The minimal setup you need to do to run a simulation is setting a domain and an array of diffusivities like this:

from diffsim.mc import MC
import numpy as np

domain = np.random.randint(0, 4, (100, 100, 100), dtype=np.uint32)
diffusivities = np.array([1e-9, 2e-10, 2e-9, 3e-9], dtype=np.float32)
mc = MC(domain=domain, diffusivities=diffusivities)
result = mc.run()

Attributes can also be updated after the MC object has been created or in the run method. This can be useful if you want to run multiple simulations with the same domain but different parameters. For example:

mc = MC(domain=domain, diffusivities=diffusivities)
result1 = mc.run(n_walkers=100000)
mc.diffusivities = np.array([1e-9, 2e-10, 2e-9, 3e-9], dtype=np.float32)
result2 = mc.run(n_walkers=200000)

Results are returned in a MCRun object. See here for documentation on that.

!!! warning The CUDA kernel that gets launched if device=gpu is not interruptable from the CPU, which means that you cannot stop it with Ctrl+C. If you want to interrupt the simulation, you need to kill the process. The easiest way to do this is to open nvtop, go to your process, press F9 (fn+F9 on Mac), choose 9 SIGKILL and press Enter.

!!! note If you want to run very large simulations with many walkers and a larger domain, you can set mc.n_runs_split to a value > 1. This will split the runs into n parts and run them sequentially. It requires less memory and can therefore be faster. This way of running a simulation will combine the results in a reconstructed array. The rest of the information will be lost (i.e. signals, xts, xis). It only changes has an effect if you run on the GPU.

Reconstruction

As each walker holds their individual signal at the end of the simulation, this information needs to be combined into an image domain again with a reconstruction. The reconstruction is done by summing up the signals of all walkers based on their final position in the domain.

Device management

By default MC runs on the GPU if available. Running on the GPU is about 50 times faster. diffsim will automatically chose a GPU for you based on availability, otherwise you can specify the GPU yourself either bei changing mc_settings.gpu_id before creating the MC object or by setting mc.gpu_id = 2. For more details on the mc_settings see here. By default MC will leave results and input data on the GPU after a simulation to avoid unnecessary data transfer. You can move all relevant data around with the cpu() and gpu() methods.

result = mc.run()
result.cpu()

Default values

If not specified otherwise, MC will use the mc_settings parameter for device management. See default values of that here.

The MC class itself has a few default values set, which will influence the simulation. See below for those default values: ::: diffsim.mc.MC options: show_root_heading: true members:

n_walkers
dt
voxel_size
do_reconstruction
reconstruct_downsample
save_intermediate_steps
save_every_x_steps
fixed_xi
seed
save_results

For a full list of available attributes have a look at the source code.

Other ways of configuring the simulation

You can also create the MC class from a yaml file containing the configuration. This can be useful if you want to run the simulation with different parameters multiple times. The yaml file should look like this:

n_walkers: 100000
dt: 1e-6
voxel_size: 1e-6

You can then create the MC object with the from_yaml method: ::: diffsim.mc.MC options: show_root_heading: true members:

from_yaml

If you want to load tissue parameters for a specific simulation, you can load a json file with the tissue parameters. The json file should look like this:

{
    "0": {
        "name": "background",
        "diffusivity": 1e-09,
        "permeability": {
            "0": 0,
            "1": 0,
            "2": 0,
        },
        "T2": 0.17
    },
    "1": {
        "name": "epithelium",
        "diffusivity": 1e-09,
        "permeability": {
            "0": 0,
            "1": 0.0011,
            "2": 0.5,
        },
        "T2": 0.1
    },
    "2": {
        "name": "lumen",
        "diffusivity": 2e-09,
        "permeability": {
            "0": 0,
            "1": 0.5,
            "2": 0.0011,
        },
        "T2": 0.14
    }
}

You can give the path to the json file like this:

mc = MC(path_to_domain_metadata="path/to/json")

Attributes and their meaning​

How to run a simulation​

Reconstruction​

Device management​

Default values​

Other ways of configuring the simulation​