Flux balance analysis (FBA) based bilevel optimisation has been a great success in redesigning metabolic networks for biochemical overproduction. To date, many computational approaches have been developed to solve the resulting bilevel optimisation problems. However, most of them are of limited use due to biased optimality principle, poor scalability with the size of metabolic networks, potential numeric issues, or low quantity of design solutions in a single run. In this work, we have employed a network interdiction model free of growth optimality assumptions, a special case of bilevel optimisation, for computational strain design and have developed a hybrid Benders algorithm (HBA) that deals with complicating binary variables in the model, thereby achieving high efficiency without numeric issues in search of best design strategies. More importantly, HBA can list solutions that meet users' production requirements during the search, making it possible to obtain numerous design strategies at a small runtime overhead (typically ∼1 hour). * math4neu@gmail.com However, the bilevel optimisation based framework in literature has numerous limitations. The first one is the intensive computational cost in search of optimal solutions. Bilevel optimisation is often reformulated into a mixed-integer linear programming (MILP) so as to be solved by exact MILP solvers. It can take up to a week to solve an MILP resulting from a standard GSMM [Feist et al., 2010]. Many practical strategies, such as model reduction and refinement of candidate knockout set [Feist et al., 2010], have been used to reduce the computational time but may miss the best design strategies due to reduced search space. GDBB introduced a truncated branch and bound to speed up the search process. GDLS used local search with multiple search paths to reduce the search space for each local MILP [Lun et al., 2009]. While finding optimal solutions are computationally costly for exact solvers, other studies resorts to inexact methods, such as genetic algorithms [Patil et al., 2005;Rocha et al., 2010] and swarm intelligence [Choon et al., 2015]. These methods, however, still scale poorly with the size of GSMM and are specially ineffective when a large number of genetic manipulations are allowed for high target production.In company with intensive computations, the resulting MILP often has weak LP relaxations due to disjunctive big-M constraints [Codato and Fischetti, 2006], another limitation of the current bilevel optimisation based methods. Big-M formulation can very easily cause numeric issues, particularly in genome-scale metabolic models where stoichiometric coefficients often vary many orders of magnitude [Sun et al., 2013]. As a result, optimal strain design solutions returned from global search solvers such as Gurobi and CPLEX may turn out to be in silico infeasible. Model reformulation may alleviate numeric issues and potentially reduce computational costs. However, a proper model reformulation is often time-consuming and laborious as extra care has to be ta...