code: def kurtosis_element(
      D_comps, frac, ind_i, ind_j, ind_k, ind_l, DT=None, MD=None
  ):
      r"""Computes the diffusion kurtosis tensor element (with indexes i, j, k
      and l) based on the individual diffusion tensor components of a
      multicompartmental model.
  
      Parameters
      ----------
      D_comps : (K,3,3) ndarray
          Diffusion tensors for all K individual compartment of the
          multicompartmental model.
      frac : [float]
          Percentage of the contribution of each tensor. The sum of fractions
          should be equal to 100%.
      ind_i : int
          Element's index i (0 for x, 1 for y, 2 for z)
      ind_j : int
          Element's index j (0 for x, 1 for y, 2 for z)
      ind_k : int
          Element's index k (0 for x, 1 for y, 2 for z)
      ind_l: int
          Elements index l (0 for x, 1 for y, 2 for z)
      DT : (3,3) ndarray (optional)
          Voxel's global diffusion tensor.
      MD : float (optional)
          Voxel's global mean diffusivity.
  
      Returns
      -------
      wijkl : float
              kurtosis tensor element of index i, j, k, l
  
      Notes
      -----
      wijkl is calculated using equation 8 given in [1]_
  
      References
      ----------
      .. [1] R. Neto Henriques et al., "Exploring the 3D geometry of the
             diffusion kurtosis tensor - Impact on the development of robust
             tractography procedures and novel biomarkers", NeuroImage (2015)
             111, 85-99.
  
      """
      if DT is None:
          DT = np.zeros((3, 3))
          for i in range(len(frac)):
              DT = DT + frac[i] * D_comps[i]
  
      if MD is None:
          MD = (DT[0][0] + DT[1][1] + DT[2][2]) / 3
  
      wijkl = 0
  
      for f in range(len(frac)):
          wijkl = wijkl + frac[f] * (
              D_comps[f][ind_i][ind_j] * D_comps[f][ind_k][ind_l]
              + D_comps[f][ind_i][ind_k] * D_comps[f][ind_j][ind_l]
              + D_comps[f][ind_i][ind_l] * D_comps[f][ind_j][ind_k]
          )
  
      wijkl = (
          wijkl
          - DT[ind_i][ind_j] * DT[ind_k][ind_l]
          - DT[ind_i][ind_k] * DT[ind_j][ind_l]
          - DT[ind_i][ind_l] * DT[ind_j][ind_k]
      ) / (MD**2)
  
      return wijkl
  
imports: import numpy
  from numpy import dot
  from dipy.core.geometry import sphere2cart
  from dipy.core.geometry import vec2vec_rotmat
  from dipy.core.gradients import GradientTable
  from dipy.reconst.utils import dki_design_matrix
  from scipy.special import jn
label_desc: Use image node_id 13 for calling a global function with example usage: kurtosis_element(D_comps, frac, ind_i, ind_j, ind_k, ind_l, DT, MD) and returns: wijkl