The main emphasis of this paper is the investigation of propellant-minimizing interplanetary low-thrust trajectories modeled through differential inclusion constraints in terms of the modified equinoctial orbital elements. An additional component is the use of analytic gradients of the transcribed equations of motion. An Earth-to-Mars test case is presented to demonstrate that the implementation of the higher-order derivative approximations and the new differential constraints is valid. Several missions including both nuclear electric and solar electric propulsion systems are described. The results show good agreement between this new method and solutions obtained with industry-standard software.