Summary: We consider a model for gene regulatory networks that is a modification of S. A. Kauffman’s [J. Theor. Biol. 22, 437–467 (1969)] random Boolean networks. There are three parameters: \(n=\) the number of nodes, \(r=\) the number of inputs to each node, and \(p=\) the expected fraction of 1’s in the Boolean functions at each node. Following a standard practice in the physics literature, we use a threshold contact process on a random graph on \(n\) nodes, in which each node has in degree \(r\), to approximate its dynamics. We show that if \(r\geq 3\) and \(r\cdot 2p(1- p)> 1\), then the threshold contact process persists for a long time, which correspond to chaotic behavior of the Boolean network. Unfortunately, we are only able to prove the persistence time is \(\geq\exp(cn^{b(p)})\) with \(b(p)> 0\) when \(r\cdot 2p(1- p)> 1\), and \(b(p)= 1\) when \((r- 1)\cdot 2p(1- p)> 1\).