My, very limited, understanding of the problem is that the Geodic system is a mathematical model where as the MSL are local observations which inherently take into account numerous factors that are not in the model. Therefore any correlation/transformation would need to be the form of localized correction. What I'm trying to say is I don't think you will find a formula that will convert from on to the other for all points.
Follow up to my last reply. The British Ordinance survey group, http://www.ordnancesurvey.co.uk, use the Airy ellipsoid model to propagate MSL from the Newlyn Datum. They have their own transformation model (OSGN02) to transform between Airy & WSG84. This is only applicable to the Newlyn Datum. I assume other countries will have similar transforms for their datums
World-wide, various standard 'geoids' approximate MSL and, in turn, these are approximated themselves as the sum of spherical harmonics with coefficients determined by actual MSL measurements. So in this way it is possible to actually calculate the localised height differences between the WGS84 ellipsoid and the approximating MSL geoid. The calculation is pretty involved so if that's what you want to do, then try using something like GDAL (http://www.gdal.org/) or other open sources (I found http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html). The theory can be found in a number of places (http://earth-info.nga.mil/GandG/publications/tr8350.2/tr8350.2-a/Chapter%206.pdf). Hope it's what you need.