How to interpret the GROMACS MD Simulation RMSD Graph? Hello brilliant researchers!
After virtual screening, I performed MD Simulation via GROMACS to check the stability of the protein-top hit complex by following this tutorial http://www.mdtutorials.com/gmx/complex/09_analysis.html

After that, I facing some difficulties in interpreting the resulted RMSD graph as it is quite different from most of the published RMSD graph. 

Based on the graph attached, how to interpret it? Is AChE-dexketoprofen complex (Red) more stable (or conserved) than BuChE-dexketoprofen complex (Black)?
The opposite-ish.

Your graph is showing only the ligand, but RMSD is generally used for structures as a whole —so I'll mention about whole protein for future readers.

RMSD is a measure of how much does a structure deviate in position from a reference. Your ligand in the protein in red "wiggles" more than the one in black.

Structure RSMD

If this were RMSD of the protein (and not the index group for the ligand), it would not really mean that the black protein is more stable. Stability is to do with force-field energetics not position. Some protein have two state positions, say MBP is a flexing dumbbell. This would be seen as a jump in the RSMD graph, similar to that seen for the black trace.

However, the command you used was based on gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg from the page you linked. Key here is the index -n index.ndx, which limits it to the ligand. Before this happened gmx trjconv aligned the protein in the trajectory, so that common mistake was not done.

Ligand RSMD

In terms of ligand RMSD, lower deviation/less wobble correlates with more stable binding, but energetics (gmx energy) will tell you how strongly does it bind.

The black trace has a jump, so your ligand in the black trace (which wobbles less) has an alternative binding mode. This is why you are looking at RSMD.
To explore this further you can:


look at the trajectory
in PyMOL using the command split_states twice, once for the first and then the last state, followed by align of the two
If it were the whole protein, looking at the per bond matrix -bm option in Gromacs rms and spotting the change.


3Å??

I strongly suggest checking the trajectory because both have a deviation of 3Å, which for a ligand is not good.


Did the energy minimisation converge properly?
If it is because the proteins themselves wobble substantially, revisit gmx trjconv with a focus on the active site minus the ligand.
Is the ligand binding but leaving one substituent group to wobble?


Energetics

As with docking, the "score" is free energy difference (in kcal/mol). This is not the same as k_off and k_on and their ratio K_d (Dissociation constant). For those you need to do steered MD or have a way longer simulation (not recommended). The reason being is that the ligand has to make its way into the active site. However, you have a secondary binding mode in black, so have two free energy scores, this means that there is a stronger k_on.