docs: def run(self):
          """Run flux variability command"""
  
          # Load compound information
          def compound_name(id):
              if id not in self._model.compounds:
                  return id
              return self._model.compounds[id].properties.get('name', id)
  
          reaction = self._get_objective()
          if not self._mm.has_reaction(reaction):
              self.fail(
                  'Specified reaction is not in model: {}'.format(reaction))
  
          loop_removal = self._get_loop_removal_option()
          enable_tfba = loop_removal == 'tfba'
          if enable_tfba:
              solver = self._get_solver(integer=True)
          else:
              solver = self._get_solver()
  
          start_time = time.time()
  
          try:
              fba_fluxes = dict(fluxanalysis.flux_balance(
                  self._mm, reaction, tfba=False, solver=solver))
          except fluxanalysis.FluxBalanceError as e:
              self.report_flux_balance_error(e)
  
          threshold = self._args.threshold
          if threshold.relative:
              threshold.reference = fba_fluxes[reaction]
  
          logger.info('Setting objective threshold to {}'.format(
              threshold))
  
          handler_args = (
              self._mm, solver, enable_tfba, float(threshold), reaction)
          executor = self._create_executor(
              FVATaskHandler, handler_args, cpus_per_worker=2)
  
          def iter_results():
              results = {}
              with executor:
                  for (reaction_id, direction), value in executor.imap_unordered(
                          product(self._mm.reactions, (1, -1)), 16):
                      if reaction_id not in results:
                          results[reaction_id] = value
                          continue
  
                      other_value = results[reaction_id]
                      if direction == -1:
                          bounds = value, other_value
                      else:
                          bounds = other_value, value
  
                      yield reaction_id, bounds
  
              executor.join()
  
          for reaction_id, (lower, upper) in iter_results():
              rx = self._mm.get_reaction(reaction_id)
              rxt = rx.translated_compounds(compound_name)
              print('{}\t{}\t{}\t{}'.format(reaction_id, lower, upper, rxt))
  
          logger.info('Solving took {:.2f} seconds'.format(
              time.time() - start_time))