The complexity of gene regulatory networks in terms of both large-scale description as well as nonlinear models is often an obstacle for analysis purposes. Therefore, the development of eective model reduction techniques is of paramount importance in the eld of systems biology. In this paper, we apply Carleman biliniarization for model reduction for gene regulatory networks based only on gramians computations. The method is based on the bilinear representation of weakly nonlinear systems and Taylor's series expansion. Thus, we obtain a computational simple solution and identify parameters that are relevant to the behavior of the system. The theoretical results are elucidated in an illustrative example and thus shown how they can be applied to reverse engineering design.