Simulation Module

Inputs and Launchers

Parameters Class

Simulation Class

Setup Functions

spin_init_positions._find_spin_locations(self)[source]

Helper function to find initial spin locations

spin_init_positions._find_spin_locations_kernel(resident_fiber_indxs_cuda: DeviceNDArray, resident_cell_indxs_cuda: DeviceNDArray, fiber_centers_cuda: DeviceNDArray, fiber_directions_cuda: DeviceNDArray, fiber_radii_cuda: DeviceNDArray, cell_centers_cuda: DeviceNDArray, cell_radii_cuda: DeviceNDArray, spin_positions_cuda: DeviceNDArray, theta_cuda, curvature_params) → None[source]

Locate spins within resident microstructural elements

Parameters:

  • resident_fiber_indxs_cuda (CUDA Device Array) – Array to write resident fiber indices to

  • resident_cell_indxs_cuda (CUDA Device Array) – Array to write resident cell indices to

  • fiber_centers_cuda (CUDA Device Array) – Coordinates for the centers of each fiber

  • fiber_directions_cuda (CUDA Device Array) – Direction of each fiber

  • fiber_radii_cuda (CUDA Device Array) – Radii of each fiber

  • cell_centers_cuda (CUDA Device Array) – Coordinates for the centers of each cell

  • cell_radii_cuda (CUDA Device Array) – Radii of each cell

  • spin_positions_cuda (CUDA Device Array) – Initial spin positions

Class Objects

Cells

class objects.cell(cell_center, cell_radius: float, cell_diffusivity: float)[source]
__init__(cell_center, cell_radius: float, cell_diffusivity: float) → None[source]

Cell information

Parameters:

  • cell_center (np.ndarray) – (x,y,z) coordinates of the cell center

  • cell_radius (float) – Intrinsic fiber diffusivity

  • cell_diffusivity (float) – Fiber radius

_get_center()[source]

Returns an array containing coordinates for the center of the specified cell

Returns:

An array containing the cell center coordinates

Return type:

numpy.ndarray

_get_diffusivity()[source]

Returns the diffusivity within the specified cell in units of \({\mathrm{μm}^2}\, \mathrm{ms}^{-1}\)

Returns:

The diffusivity within the cell

Return type:

float

_get_radius()[source]

Returns the radius of the specified cell, in units of \({\mathrm{μm}}\)

Returns:

The cell radius

Return type:

float

_set_center(cell_center: ndarray)[source]

Records the center coordinates for the specified cell

Parameters:

cell_center (np.ndarray) – An array containing the cell center coordinates

_set_diffusivity(D0: float)[source]

Records the diffusivity within the specified cell in units of \({\mathrm{μm}^2}\, \mathrm{ms}^{-1}\)

Parameters:

D0 (float) – The diffusivity within the cell, assumed to be equal to the user-specified free water diffusivity

_set_radius(radius: float)[source]

Records the radius of the specified cell, in units of \({\mathrm{μm}}\)

Parameters:

radius (float) – The cell radius

Fibers

class objects.fiber(center: ndarray, direction: ndarray, bundle: int, diffusivity: float, radius: float, kappa: float, L: float, A: float, P: float)[source]

Class object for fiber attributes

__init__(center: ndarray, direction: ndarray, bundle: int, diffusivity: float, radius: float, kappa: float, L: float, A: float, P: float) → None[source]

Fiber information and parameters

Parameters:

  • center (np.ndarray) – (x,y,z) coordinates of the fiber center

  • direction (np.ndarray) – Constituent bundle index

  • bundle (int) – Unit vector pointed along the fiber direction

  • diffusivity (float) – Intrinsic fiber diffusivity, square micrometers per millisecond

  • radius (float) – Fiber radius, in micrometers

_get_bundle()[source]

Returns the fiber bundle index for the specified fiber

Returns:

The bundle/type index for the specified fiber

Return type:

int

_get_center()[source]

Returns an array containing coordinates for the center of the specified fiber

Returns:

An array of fiber center coordinates

Return type:

numpy.ndarray

_get_diffusivity()[source]

Returns the intra-axonal diffusivity for the specified fiber, in units of \({\mathrm{μm}^2}\, \mathrm{ms}^{-1}\)

Returns:

The user-specified diffusivity within the specified fiber

Return type:

float

_get_direction()[source]

Returns the orientation vector for the specified fiber

Returns:

The orientation vector/rotation matrix for the specified fiber

Return type:

numpy.ndarray

_get_radius()[source]

Returns the radius of the specified fiber, in units of \({\mathrm{μm}}\)

Returns:

The user-specified radius of the specified fiber

Return type:

float

Spins

class objects.spin(spin_position_t1m: ndarray)[source]
__init__(spin_position_t1m: ndarray) → None[source]

Spin information

Parameters:

  • spin_position_t1m (np.ndarray) – Initial spin position

  • spin_position_t1m – Final spin position

  • in_fiber_index (int) – Index of resident fiber (if the spin resides in a fiber)

  • fiber_bundle (int) – Index of resident fiber bundle (if the spin resides in a fiber)

  • in_cell_index (int) – Index of the resident cell (if the spin resides in a cell)

  • in_water_index (int) – Spin index if in water

_get_bundle_index()[source]

Returns the fiber bundle index for a spin residing in a fiber

Returns:

Index of the resident fiber bundle (if the spin resides in a fiber)

Return type:

int

_get_cell_index()[source]

Returns the index of the spin if it resides in a cell

Returns:

Index of the resident cell (if the spin resides in a cell)

Return type:

int

_get_fiber_index()[source]

Returns the index of the fiber in which the spin resides.

Returns:

Index of the resident fiber bundle (if the spin resides in a fiber)

Return type:

int

_get_position_t1m()[source]

Returns the initial position of the spin

Returns:

Initial spin position

Return type:

int

_get_position_t2p()[source]

Returns the final position of the spin

Returns:

Final spin position

Return type:

int

_get_water_index()[source]

Returns the water index of the specified spin

Returns:

Spin index if in water

Return type:

int

_set_cell_index(index: int)[source]

Records the index of the spin if it resides in a cell

Parameters:

index (int) – Index of the resident cell (if the spin resides in a cell)

_set_fiber_bundle(index: int)[source]

Records the fiber bundle for a spin residing in a fiber

Parameters:

index (int) – Index of the resident fiber bundle (if the spin resides in a fiber)

_set_fiber_index(index: int)[source]

Records the index of the fiber in which the spin resides.

Parameters:

index (np.ndarray) – Index of resident fiber (if the spin resides in a fiber)

_set_position_t2p(position: ndarray)[source]

Records the position of the spin after diffusion time \(\Delta\).

Parameters:

position (np.ndarray) – Final spin position

_set_water_index(index: int)[source]

Records the index of the spin if it resides in water

Parameters:

index (int) – Spin index if in water

Diffusion Physics

Wrapper

diffusion._caclulate_volumes(spins)[source]

Calculates empirical volume fraction for each simulated compartment (i.e., fibers, cells, and water)

Parameters:

spins (list) – A one-dimensional list (length = n_walkers) containing each instance of objects.spin(). Each entry corresponds to one spin from the spin ensemble.

diffusion._diffusion_context_manager(random_states, spin_positions, spin_in_fiber_at_index, fiber_centers, fiber_step, fiber_radii, fiber_directions, spin_in_cell_at_index, cell_centers, cell_step, cell_radii, water_step, theta_cuda, void, curvature_params)[source]

Helper function to segment each spin into the relevant physics module for its resident compartment

Parameters:

  • random_states – Randomized states

  • spin_positions – Initial spin positions

  • spin_in_fiber_at_index – Spin indices for each fiber

  • fiber_centers – Geometric centers of each fiber

  • fiber_step – Length of diffusion step to take for each fiber type, in units of micrometers

  • fiber_radii – Radius of each fiber, in units of micrometers

  • fiber_directions – Relative fiber rotations

  • spin_in_cell_at_index – Spin indices for each cell

  • cell_centers – Geometric centers for each cell

  • cell_step – Length of diffusion step to take for each cell type, in units of micrometers

  • cell_radii – Radius of each cell, in units of micrometers

  • water_step – Length of diffusion step to take in free water, in units of micrometers

  • void – Boolean argument, True if fiber configuration is void, False otherwise

diffusion._simulate_diffusion(self) → None[source]

Iterates over the range \(t \in [0, \Delta ]\) with a step size of \(\dd{t}\).

Parameters:

  • self (class object) – the dmri_simulation object

  • spins (list) – A one-dimensional list (length = n_walkers) containing each instance of objects.spin(). Each entry corresponds to one spin from the spin ensemble.

  • cells (list) – A one-dimensional list (length = n_cells) containing each instance of objects.cell(). Each entry corresponds to one cell within the simulated imaging voxel.

  • fibers (list) – A list (size = n_fibers\(\times\)n_fibers) containing each instance of objects.fiber(). Each entry corresponds to one fiber within the simulated imaging voxel.

  • Delta (float) – The diffusion time supplied by the user in units of milliseconds.

  • dt (float) – The user-supplied duration for each time step in units of milliseconds. Assumed to be equal to \(\delta\) under the narrow-pulse approximation.

  • water_diffusivity (float) – The user-supplied diffusivity for free water, in units of \({\mathrm{μm}^2}\, \mathrm{ms}^{-1}\).

Returns:

Updated spin trajectories within each instance of the spin object in the spin list

Note

At each iteration, the updated spin position is written to spin_positions_cuda

Diffusion in Cells

walk_in_cell._diffusion_in_cell(i, random_states, cell_center, cell_radii, cell_step, fiber_centers, fiber_radii, fiber_directions, spin_positions, thetas, void, A)[source]

Simulated Brownian motion of a spin confined to within in a cell, implemented via random walk with rejection sampling for proposed steps beyond the cell membrane. Note that this implementation assumes zero exchange between compartments and is therefore only pysically-accurate for \(\Delta < {\tau_{i}}\) [1].

Parameters:

  • i (int) – Absolute index of the current thread within the block grid

  • random_states (numba.cuda.cudadrv.devicearray.DeviceNDArray) – xoroshiro128p random states

  • cell_center (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Coordinates of the cell centers

  • cell_radii (float) – Cell radius, in units of \({\mathrm{μm}}\)

  • cell_step (float) – Distance travelled by resident spins for each time step \(\dd{t}\)

  • fiber_centers (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Coordinates of the fiber centers

  • fiber_radii (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Radii of each fiber type

  • fiber_directions (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Orientation of each fiber type

  • spin_positions (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Array containing the updated spin positions

  • void (bool) – Logical condition that is True if fiber_configuration = Void and False otherwise

Shapes
random_states:

(n_walkers,) where n_walkers is an input parameter denoting the number of spins in the ensemble

cell_center:

(3,)

fiber_centers:

(n_fibers x n_fibers, 3) where n_fibers is computed in a manner such that the fibers occupy the supplied fiber fraction of the imaging voxel

fiber_radii:

(n_fibers x n_fibers, ) where n_fibers is computed in a manner such that the fibers occupy the supplied fiber fraction of the imaging voxel

fiber_directions:

(n_fibers x n_fibers, 3) where n_fibers is computed in a manner such that the fibers occupy the supplied fiber fraction of the imaging voxel

spin_positions:

(n_walkers, 3) where n_walkers is an input parameter denoting the number of spins in the ensemble

Diffusion in Fibers

walk_in_fiber._diffusion_in_fiber(i, random_states, fiber_center, fiber_radius, fiber_direction, fiber_step, theta, spin_positions, curvature_params)[source]

Simulated Brownian motion of a spin confined to within in a fiber, implemented via random walk with rejection sampling for proposed steps beyond the fiber membrane. Note that this implementation assumes zero exchange between compartments and is therefore only physically-accurate for \(\Delta < {\tau_{i}}\) [1].

Parameters:

  • i (int) – Absolute index of the current thread within the block grid

  • random_states (numba.cuda.cudadrv.devicearray.DeviceNDArray) – xoroshiro128p random states

  • fiber_center (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Coordinates of the center of specified fiber

  • fiber_radius (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Radius of the specified fiber type

  • fiber_direction (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Orientation of the specified fiber type

  • fiber_step (float) – Distance travelled by resident spins for each time step \(\dd{t}\)

  • spin_positions (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Array containing the updated spin positions

Shapes
random_states:

(n_walkers,) where n_walkers is an input parameter denoting the number of spins in the ensemble

cell_center:

(3,)

fiber_centers:

(n_fibers x n_fibers, 3) where n_fibers is computed in a manner such that the fibers occupy the supplied fiber fraction of the imaging voxel

fiber_radii:

(n_fibers x n_fibers, ) where n_fibers is computed in a manner such that the fibers occupy the supplied fiber fraction of the imaging voxel

fiber_directions:

(n_fibers x n_fibers, 3) where n_fibers is computed in a manner such that the fibers occupy the supplied fiber fraction of the imaging voxel

spin_positions:

(n_walkers, 3) where n_walkers is an input parameter denoting the number of spins in the ensemble

Diffusion in Water

walk_in_water._diffusion_in_water(i, random_states, fiber_centers, fiber_directions, fiber_radii, cell_centers, cell_radii, spin_positions, step, thetas, curvature_params)[source]

Simulated Brownian motion of a spin confined to the extra-cellular/axonal water, implemented via random walk with rejection sampling for proposed steps into cells or fibers. Note that this implementation assumes zero exchange between compartments and is therefore only pysically-accurate for \(\Delta < {\tau_{i}( ext{cells})}\) and \(\Delta < {\tau_{i}( ext{fibers})}\) [1].

Parameters:

  • i (int) – Absolute index of the current thread within the block grid

  • random_states (numba.cuda.cudadrv.devicearray.DeviceNDArray) – xoroshiro128p random states

  • fiber_centers (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Coordinates of the fiber centers

  • fiber_directions (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Orientation of each fiber type

  • fiber_radii (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Radii of each fiber type

  • cell_centers (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Coordinates of the cell centers

  • cell_radii (float) – Cell radius, in units of \({\mathrm{μm}}\)

  • spin_positions (numba.cuda.cudadrv.devicearray.DeviceNDArray) – Array containing the updated spin positions

  • step (float) – Distance travelled by resident spins for each time step \(\dd{t}\)

Shapes
random_states:

(n_walkers,) where n_walkers is an input parameter denoting the number of spins in the ensemble

cell_center:

(3,)

fiber_centers:

(n_fibers x n_fibers, 3) where n_fibers is computed in a manner such that the fibers occupy the supplied fiber fraction of the imaging voxel

fiber_radii:

(n_fibers x n_fibers, ) where n_fibers is computed in a manner such that the fibers occupy the supplied fiber fraction of the imaging voxel

fiber_directions:

(n_fibers x n_fibers, 3) where n_fibers is computed in a manner such that the fibers occupy the supplied fiber fraction of the imaging voxel

spin_positions:

(n_walkers, 3) where n_walkers is an input parameter denoting the number of spins in the ensemble

Save Outputs

save._add_noise(signal, snr)[source]

Add Gaussian noise to the forward simulated signal [2]

Parameters:

  • signal (np.ndarray) – The forward simulated signal

  • snr (float) – The desired signal to noise ratio of the b0 image

Returns:

The noise-added signal

Return type:

np.ndarray

save._generate_signals_and_trajectories(self)[source]

Helper function to organize and store compartment specific and combined trajectories and their incident signals

Returns:

signals with associated labels (signals_dict), trajectories with associated labels (trajectories_dict)

Return type:

dictionaries

save._save_data(self)[source]

Helper function that saves signals and trajectories to the current directory.

save._signal(phase: ndarray, SNR: int | None = None) → ndarray[source]

Calculates the PGSE signal from the forward simulated spin trajectories [3]. Note that this computation is executed on the GPU using PyTorch.

Parameters:

  • spins (list) – A list of each objects.spin instance corresponding to a spin in the ensemble of random walkers

  • bvals (np.ndarray) – The supplied b-values (diffusion weighting factors)

  • bvecs (np.ndarray) – The supplied diffusion gradients

  • Delta (float) – The diffusion time, in milliseconds

  • dt (float) – The time step parameter, also equal to delta because of the narrow pulse approximation

  • SNR (float, optional) – The snr of the b0 image. If a value is not entered, the SNR of the signal is infinite. (Defaults to None)

Returns:

the forward simulated PGSE signal (signal), the initial spin positions (trajectory_t1m), and the final spin positions trajectory_t2p

Return type:

np.ndarray

Miscellaneous Functions

linalg.Ry(thetas)[source]

Calculates rotation matrices for the desired fiber orientations

Parameters:

thetas (int, tuple) – List of angles (with respect to the y-axis) for each fiber bundle

Returns:

Rotation matrices for each fiber bundle.

Return type:

np.ndarray

linalg.affine_transformation(xv_total: ndarray, yv_total: ndarray, x: float, y: float, thetas, i)[source]

Calculates and applies affine transformation to fiber grid.

Parameters:

  • xv (np.ndarray) – Fiber grid

  • x (float) – x coordinates

  • y (float) – y coordinates

  • thetas (int, tuple) – _description_

  • i (int) – bundle index

Returns:

Transformed coordinates

Return type:

np.ndarray

linalg.dL2(x, y, v, project_along)[source]

Internal linear algebra function for projecting along transformed vectors.

Parameters:

  • x (float) – x coordinate

  • y (float) – x coordinate

  • v (list) – Vector to project along

Returns:

Distance

Return type:

float

References