The standard way to compute the time variation of the orbital elements starts from Lagrange's planetary
equations. Other answers on this site which describe them are here, here, and here. These equations don't have any specific answer; they just list the relations among variables in the two-body problem. Think of them as "F=ma" for Keplerian elements: given a set of perturbing forces, they will calculate what the change in state (the time derivatives of the elements) must be, but you have to provide the forces, and the initial state, as inputs. At first glance, it looks like the third-body gravity of the sun has to be the biggest by far, but let's not start with that.
Instead, even before trying to figure out which perturbations are most important to consider, we can see possible trouble ahead. Three of the six equations, those for ascending node $\Omega$, inclination $i$, and argument of periapsis $\omega$ contain $\sqrt{1-e^2}\sin i$ in their denominators, so as eccentricity $e$ approaches one and inclination approaches 180 degrees, all those terms get big. They don't get super huge ( $1/\sqrt{1-e^2}$ ranges from 2.4 at $e=0.91$ to 4.1 at $e=0.97$, and is not quite 7.1 even at $e=0.99$; inclination goes faster, as 1/sin 170 = 5.8, 1/sin 175 = 11.5, and 1/sin 179 = 57.3), but it's enough to be unusually large, since that eccentricity is most definitely unusually large. The $1/\!\sin i$ term tells us something breaks at $i=\pi$: the problem is $\Omega$ is undefined there (there is no ascending node, because the orbit never ascends out of the equatorial plane), so for $i$ near $\pi$ (or zero) a small change in the orbit estimate often produces a large change in $\Omega$. What breaks at $e=1$ is that it isn't an orbit anymore, because it's an escape trajectory, so nothing repeats.
There are other parameterizations, such as the equinoctial elements, that try to avoid these and other singularities and inconveniences, at the cost of being somewhat harder to draw or describe. Sometimes the answer to "why is this element changing like crazy" is "because it's always badly behaved when X happens, so you should use something else instead", but I'm not actually going to get into them here.
What I want to draw your attention to is the curves for $\Omega$ and $\omega$, which look almost identical, except for a nearly constant offset. This behavior is expected for certain types of perturbations, but not others. That might tell us something about what kinds of perturbations are strongest, but in this case there's not much else that could have happened. In particular, consider the two Lagrange planetary equations
$$ \frac{d\Omega}{dt} = \frac{1}{nab\sin i}\frac{\partial R}{\partial i} $$
$$ \frac{d\omega}{dt} = \frac{-\cos i}{nab\sin i}\frac{\partial R}{\partial i} + \frac{b}{na^3e}\frac{\partial R}{\partial e}$$
where $a$ is the semi-major axis, $b = a\sqrt{1-e^2}$ is the semi-minor axis, $n = \sqrt{GM/a^3}$ is the mean motion, and $R$ is called "the disturbing function". How it is chosen to represent different perturbations is more complicated than I want to get into here. What I want to point out is, first, that in this orbit, since $(-\!\cos i)$ is basically one, the first term in $\dot{\omega}$ pretty much equals $\dot{\Omega}$. Second, the ratio of the coefficients of the two partials of $R$ in $\dot{\omega}$ is
$$ \frac{b}{na^3e} \div \frac{-\cos i}{nab\sin i} = \frac{nab^2 \sin i}{-na^3e \cos i} = \frac{(1-e^2)\sin i}{-e \cos i} \approx \frac{0\times 0}{1 \times 1}$$
so the contribution of $\frac{\partial R}{\partial e}$ to $\dot{\omega}$ is almost certainly ignorable when compared with that of $\frac{\partial R}{\partial i}$. This combination of $e$ and $i$ pretty much requires that $\Omega=\omega\ +$ constant. Also note that the step ups in the omegas are when $i$ gets closest to $\pi$, so $1/\!\sin i$ gets closest to dividing by zero, so the omega dots get much bigger rather suddenly but only for a short time.
There's much more that could be said, but I'll stop here. I also think some time series analysis, particularly a power spectrum of the eccentricity curve, would highlight two dominant time scales we should be able to relate to physical periods, such as the length of the Venusian year.
