The aim of the present study is to analyze the saturation regime of the energetic-ion-driven resistive interchange mode (EIC) in the LHD plasma. A set of non linear simulations are performed by the FAR3d code that uses a reduced MHD model for the thermal plasma coupled with a gyrofluid model for the energetic particles (EP) species. The hellically trapped EP component is introduced through a modification of the averaged drift velocity operator to include their precessional drift. The non linear simulation results show similar 1/1 EIC saturation phases with respect to the experimental observations, reproducing the enhancement of the n/m = 1/1 resistive interchange modes (RIC) amplitude and width as the EP β increases, the EP β threshold for the 1/1 EIC excitation, the further destabilization of the 1/1 EIC as the population of the helically trapped EP increases and the triggering of burst events. The frequency of the 1/1 EIC calculated during the burst event is 9.4 kHz and the 2/2 and 3/3 overtones are destabilized, consistent with the frequency range and the complex mode structure measured in the experiment. In addition, the simulation shows the inward propagation of the 1/1 EIC due to the non linear destabilization of the 3/4 and 2/3 EPMs, leading to the partial overlapping between resonances during the burst event. Finally, the analysis of the 1/1 EIC stabilization phase shows the excitation of the 1/1 RIC as soon as the flattening induced by the 1/1 EIC in the pressure profile vanishes, leading to the retrieval of the pressure gradient at the plasma periphery and the overcoming of the RIC stability limit.