You might have heard about recent efforts to inspect lots of "wide binaries", double stars that orbit each other at very large distances, which is one of the tasks the Gaia mission was built for, to determine if their dynamics follows Newtonian gravity or rather MOND, the modified Newtonian dynamics (Einstein theory plays no role at such weak fields).
You can learn about the latest update from this video by Dr. Betty (spoiler: Newton's just fine).
MOND is an alternative theory of gravity that was originally proposed as an alternative to dark matter to explain galactic rotation curves (which it does quite well, some argue better than dark matter). Since, it has been investigated in other weak gravity situations as well. In short, it introduces an additional scale \(a_0\) of dimension acceleration and posits that gravitational acceleration (either in Newton's law of gravity or in Newton's second law) are weakened by a factor
$$\mu(a)=\frac{a}{\sqrt{a^2+a_0^2}}$$
where a is the acceleration without the correction.
In the recent studies reported on in the video, people measure the stars' velocities and have to do statistics because they don't know about the orbital parameters and the orientation of the orbit relative to the line of sight.
That gave me an idea of what else one could try: When the law of gravity gets modified from its \(1/r^2\) form for large separations and correspondingly small gravitational accelerations, the orbits will no longer be Keppler ellipses. What happens for example if this modified dynamics would result for example in eccentricities growing or shrinking systematically? Then we might observe too many binaries with large/small eccentricities and that would be in indication of a modified gravitational law.
The only question is: What does the modification result in? A quick internet search did not reveal anything useful combining celestial mechanics and MOND, so I had to figure out myself. Inspection shows that you can put the modification into a modification of \(1/r^2\) into
$$\mu(1/r^2) \frac{\vec r}{r^3}$$
and thus into a corresponding new gravitational potential. Thus much of the usual analysis carries over: Energy and angular momentum would still be conserved and one can go into the center of mass system and work with the reduced mass of the system. And I will use units in which \(GM=1\) to simplify calculations.
The only thing that will no longer be conserved is the RungeLenzvector
$$\vec A= \vec p\times\vec L  \vec e_r.$$
\(\vec A\) points in the direction of the major semiaxis and its length equals the eccentricity of the ellipse.
Just recall that in Newton gravity, this is an additional constant of motion (which made the system \(SO(4,2)\) rather than \(SO(3)\) symmetric and is responsible for states with different \(\ell\) being degenerate in energy for the hydrogen atom), as one can easily check
$$\dot{\vec A} = \{H, \vec A\}= \dot{\vec p}\times \vec L\dot{\vec e_r}=\dots=0$$
using the equations of motion in the first term.
To test this idea I started Mathematica and used the numerical ODE solver to solve the modified equations of motion and plot the resulting orbit. I used initial data that implies a large eccentricity (so one can easily see the orientation of the ellipse) and an \(a_0\) that kicks in for about the further away half of the orbit.
Orbit of would be Runge Lenz vector \(\vec A\) 
What a disappointment! Even if it is no longer conserved it seems to move on a circle with some additional wiggles on it (Did anybody mention epicycles?). So it is only the orientation of the orbit that changes with time but there is no general trend toward smaller or larger eccentricities that one might look out for in real data.
On the other hand the eccentricity \(\\vec A\\) is not exactly conserved but wiggles a bit with the orbit but comes back to its original value after one full rotation. Can we understand that analytically?
To this end, we make use the fact that the equation of motion is only used in the first term when computing the time derivative of \(\vec A\):
$$\dot{\vec A}=\left(1\mu(1/r^2)\right) \dot{\vec e_r}.$$
\(\mu\) differs from 1 far away from the center, where the acceleration is weakest. On the other hand, since \(\vec e_r\) is a unit vector, its time derivative has to be orthogonal to it. But in the far away part of the the ellipse, \(\vec e_r\) is almost parallel to the major semi axis and thus \(\vec A\) and thus \(\dot{\vec a}\) is almost orthogonal to \(\vec A\). Furthermore, due to the reflection symmetry of the ellipse, the parts of \(\dot{\vec e_r}\) that are not orthogonal to \(\vec A\) will cancel each other on both sides and thus the wiggling around the average \(\\vec a\\) is periodic with the period of the orbit. q.e.d.
There is only a tiny net effect since the ellipse is not exactly symmetric but precesses a little bit. This can be seen when plotting \(\\vec A\\) as a function of time:

