docs: def pathway(F, A, B):
      r"""Compute the dominant reaction-pathway.
  
      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
  
      Returns
      -------
      w : list
          The dominant reaction-pathway
  
      """
      if F.nnz == 0:
          raise PathwayError('no more pathways left: Flux matrix does not contain any positive entries')
      b1, b2, F = find_bottleneck(F, A, B)
      if np.any(A == b1):
          wL = [b1, ]
      elif np.any(B == b1):
          raise PathwayError(("Roles of vertices b1 and b2 are switched."
                              "This should never happen for a correct flux network"
                              "obtained from a reversible transition meatrix."))
      else:
          wL = pathway(F, A, [b1, ])
      if np.any(B == b2):
          wR = [b2, ]
      elif np.any(A == b2):
          raise PathwayError(("Roles of vertices b1 and b2 are switched."
                              "This should never happen for a correct flux network"
                              "obtained from a reversible transition meatrix."))
      else:
          wR = pathway(F, [b2, ], B)
      return wL + wR