docs: def pathways(F, A, B, fraction=1.0, maxiter=1000, tol=1e-14):
      r"""Decompose flux network into dominant reaction paths.
  
      Parameters
      ----------
      F : (M, M) scipy.sparse matrix
          The flux network (matrix of netflux values)
      A : array_like
          The set of starting states
      B : array_like
          The set of end states
      fraction : float, optional
          Fraction of total flux to assemble in pathway decomposition
      maxiter : int, optional
          Maximum number of pathways for decomposition
      tol : float, optional
          Floating point tolerance. The iteration is terminated once the
          relative capacity of all discovered path matches the desired
          fraction within floating point tolerance
  
      Returns
      -------
      paths : list
          List of dominant reaction pathways
      capacities: list
          List of capacities corresponding to each reactions pathway in paths
  
      References
      ----------
      .. [1] P. Metzner, C. Schuette and E. Vanden-Eijnden.
          Transition Path Theory for Markov Jump Processes.
          Multiscale Model Simul 7: 1192-1219 (2009)
  
      """
      F, a, b = add_endstates(F, A, B)
      A = [a, ]
      B = [b, ]
  
      """Total flux"""
      TF = F.tocsr()[A, :].sum()
  
      """Total capacity fo all previously found reaction paths"""
      CF = 0.0
      niter = 0
  
      """List of dominant reaction pathways"""
      paths = []
      """List of corresponding capacities"""
      capacities = []
  
      while True:
          """Find dominant pathway of flux-network"""
          try:
              path = pathway(F, A, B)
          except PathwayError:
              break
          """Compute capacity of current pathway"""
          c = capacity(F, path)
          """Remove artifical end-states"""
          path = path[1:-1]
          """Append to lists"""
          paths.append(np.array(path))
          capacities.append(c)
          """Update capacity of all previously found paths"""
          CF += c
          """Remove capacity along given path from flux-network"""
          F = remove_path(F, path)
          niter += 1
          """Current flux numerically equals fraction * total flux or is
          greater equal than fraction * total flux"""
          if (abs(CF/TF - fraction) <= tol) or (CF/TF >= fraction):
              break
          if niter > maxiter:
              warnings.warn("Maximum number of iterations reached", RuntimeWarning)
              break
      return paths, capacities