The use of machine learning (ML) for modeling is on the rise. In the age of big data, this technique has shown great potential to describe complex physical phenomena in the form of models. More recently, ML has frequently been used for turbulence modeling while the use of this technique for combustion modeling is still emerging. Gene expression programming (GEP) is one class of ML that can be used as a tool for symbolic regression and thus improve existing algebraic models using high-fidelity data. Direct numerical simulation (DNS) is a powerful candidate for producing the required data for training GEP models and validation. This paper therefore presents a highly efficient DNS solver known as HiPSTAR, originally developed for simulating non-reacting flows in particular in the context of turbo-machinery. This solver has been extended to simulate reacting flows. DNSs of two turbulent premixed jet flames with different Karlovitz numbers are performed to produce the required data for training. GEP is then used to develop algebraic flame surface density models in the context of large-eddy simulation (LES). The result of this work introduces new models which show excellent performance in prediction of the flame surface density for premixed flames featuring different Karlovitz numbers.