You could do essentially what that paper (that you linked) is doing. Their process of adding tracers (points) at/near mid-ocean ridges and following them through plate splits and merges (ie, changes in plate ID) until they subduct (or hit continents) uses essentially the same idea that GPlates currently uses (in the process I mentioned in a post above) to reconstruct arbitrary points using dynamic polygons.
The main difference is they are continually generating points along ridges through time and keeping track of the ages of these points as they are reconstructed, and then using those points and their ages to create age grids. With GPlates you decide which points to give it and when, and then it reconstructs/subducts them. In both these cases working with points is easier than lines (eg, isochrons).
The next release of pyGPlates will also be able to do this (what GPlates currently does), in which case you could then … deep breath… generate isochron lines (like the first pygplates sample code you listed), with the mid-ocean ridges coming from your own ridges (instead of topological plate boundaries like usual but the general idea is similar), and tessellate those lines into points (using pyGPlates), and move them slightly off the mid-ocean ridge (ie, make two copies so one moves on one flank of MOR and the other copy moves on other flank - you don’t want them right on the MOR making their reconstruction ambiguous) then use the yet-to-be-release pyGPlates to reconstruct them through plate splits/merges (and eventually subduct them). And you also keep generating new MOR points each time step, with any point’s age at any particular time being the point’s begin time (when created at MOR) minus the current time. Then at each time step gather all points created by MORs so far, and their ages, and create a single feature with a scalar coverage (like Simon shows in this post but it’ll be age instead of crustal thickness) - a scalar coverage is basically a feature containing a bunch of points where each point has its own scalar value(s), in this case age. And make it so that feature only exists for that particular time (ie, has [begin, end] time range of [t+0.5, t-0.5]) - it’s basically a snapshot of that time since you don’t want it to appear for any other time except t. You’ll then have one such feature for each time, together representing snapshots of the ages of all the tracked points at all times. Put all those snapshot features into a single feature collection (file) and load that into GPlates. Then you can change the time slider and see the different ages of all the points (coloured by a colour palette) through time. That should give you a pretty good visual representation of the ocean crust age through time.
That’s all contingent on the new pyGPlates release though, but I could probably send you a premature version in February. The actual public release would definitely be later than Feb though, since there’s a couple other deformation features still need to be implemented.