In this study, we investigate the boundedness, persistence of positive solutions, local and global stability of the unique positive equilibrium point and rate of convergence of positive solutions of the following difference equations systems of exponential forms:

$${\mathrm{{\rm Y}}}_{n+1}=\frac{{\mathrm{\Gamma }}_1+{\delta }_1{e}^{-{\mathrm{\Psi }}_{n-1}}}{{\mathrm{\Theta }}_1+{\mathrm{\Psi }}_n},{\mathrm{\Psi }}_{n+1}=\frac{{\mathrm{\Gamma }}_2+{\delta }_2{e}^{-{\mathrm{\Omega }}_{n-1}}}{{\mathrm{\Theta }}_2+{\mathrm{\Omega }}_n},{\mathrm{\Omega }}_{n+1}=\frac{{\mathrm{\Gamma }}_3+{\delta }_3{e}^{-{\mathrm{{\rm Y}}}_{n-1}}}{{\mathrm{\Theta }}_3+{\mathrm{{\rm Y}}}_n},$$

$${\mathrm{{\rm Y}}}_{n+1}=\frac{{\mathrm{\Gamma }}_1+{\delta }_1{e}^{-{\mathrm{\Psi }}_{n-1}}}{{\mathrm{\Theta }}_1+{\mathrm{{\rm Y}}}_n},{\mathrm{\Psi }}_{n+1}=\frac{{\mathrm{\Gamma }}_2+{\delta }_2{e}^{-{\mathrm{\Omega }}_{n-1}}}{{\mathrm{\Theta }}_2+{\mathrm{\Psi }}_n},{\mathrm{\Omega }}_{n+1}=\frac{{\mathrm{\Gamma }}_3+{\delta }_3{e}^{-{\mathrm{{\rm Y}}}_{n-1}}}{{\mathrm{\Theta }}_3+{\mathrm{\Omega }}_n},$$

for $n\in {\mathbb{N}}_0$, where the initial conditions ${\mathrm{{\rm Y}}}_{-j}$, ${\mathrm{\Psi }}_{-j}$, ${\mathrm{\Omega }}_{-j}$, for $j\in \{0,1\}$ and the parameters ${\mathrm{\Gamma }}_i$, ${\delta }_i$, ${\mathrm{\Theta }}_i$ for $i\in \{1,2,3\}$ are positive constants.