The lac operen in E. coli is one of the earliest examples of an inducible system of genes being under both positive and negative control that is capable of showing bistability. In this paper, we present a methodology to infer synthetic threshold Boolean regulatory networks of a reduced model of the lac operon using evolutionary computation. The formulation consists in a vector representation of the solutions (networks) and a fitness function specially designed to correctly simulate the bistability through the models' fixed points. We compared the effectiveness and efficiency (runtime) of the proposed approach using three evolutionary computation techniques: differential evolution, genetic algorithms, and particle swarm optimization. The results showed that the three algorithms are capable of finding solutions, being differential evolution the most effective, whereas genetic algorithms was the least effective and efficient in terms of runtime. Particle swarm optimization obtained a good trade-off between effectiveness versus efficiency. One of the inferred solutions was analyzed showing some interesting biological insights, as well as correctly being able to model bistability without any spurious attractors. Overall, the proposed formulation was effective to infer bistable lac operon models under the threshold Boolean network paradigm.