First of all, here is a program which simulates the evolution of the G-matrix over multiple generations, it's a few years old (they seem to have stopped developing it) and I've only played with it briefly. This could solve how to model the evolution of the G-matrix.
Fisher's fundamental theorem is a great place to start off with the theory of this:
The rate of increase in fitness of any organism at any time is equal to its genetic variance in fitness at that time.
What this means (as I am sure you realise but I will put in so the answer can help others too) is that evolution by selection depends not only on the strength and form of selection, but the genetic variation underlying the selected trait. This is captured in the breeders equation $\Delta \bar z = G \beta$ where $\Delta \bar z$ is the response (in a multivariate space), $G$ is a matrix of the genetic variance within traits covariance between traits (different traits, the same traits in either sex, or the same traits in multiple environments) and $\beta$ is the vector of selection gradients on all of those contexts. It seems reasonable to expect that given enough time selection will erode variation because standing genetic variance is a finite resource and selection removes polymorphism while not adding any new variants (that's for mutation & migration).
This paper discusses the effects of both selection and drift on the G-matrix, they also deal with some modelling to support their results. More specific to simulation, this article by Arnold et al (a number of the big players in the G-matrix circle). It reviews "empirical, analytical, and simulation studies of the G-matrix with a focus on its stability and evolution." It would be a really good read for you on this subject.
This paragraph captures the essence of your question:
Focusing on a longer time scale, we find that the G-matrix evolves in expected ways to the AL [adaptive landscape] and the pattern of mutation. In the absence of correlational selection (rω = 0) and mutational correlation (rμ = 0), the average G-ellipse is nearly circular, although the ellipse fluctuates wildly about this average (first row in Fig. 6). At the opposite extreme, when the leading eigenvectors of the AL and the M-matrix [mutation matrix] are both inclined at an angel of 45°, the leading eigenvector of G is pitched at the same angle (last row in Fig. 6). Between these two extremes, G tends to evolve to a shape and orientation that represents an intermediate compromise between the AL and M. In other words, the simulation results confirm our intuition and theoretical expectations (Lande 1980b) that G should evolve toward alignment with the AL and M.
Simulation studies of the G-matrix:
- Bürger R, Wagner GP, Stettinger J. How much heritable variation can be maintained in finite populations by mutation-selection balance. Evolution. 1989;43:1748–1766.
- Multivariate mutation-selection balance with constrained pleiotropic effects.
Genetics. 1989 May; 122(1):223-34.
- On the distribution of the mean and variance of a quantitative trait under mutation-selection-drift balance.
Bürger R, Lande R
Genetics. 1994 Nov; 138(3):901-12.
- Predicting long-term response to selection.
Genet Res. 2000 Feb; 75(1):83-94.
The last article there by Reeve is probably the best paper for your question because it describes in some detail the simulation model they use and how it is all set up. Briefly they simulate a population of 4000 diploid individuals with three genetically correlated traits with separate but identical sexes, random mating, and discrete generations. 20000 generations are simulated to allow mutation-selection-drift equilibrium (almost like a burn-in time on an MCMC chain). They then shift the optimum for one trait by 10 standard deviations and simulate 1500 generations in five replicates. There are 100 unlinked loci underlying the traits, with 50 loci affecting each trait (randomly assigned) meaning there is likely some degree of genetic correlation though it is not perfectly correlated (for some reading on genetic correlations and the evolution of differences you can read Bonduriansky & Rowe 2005, Poissant et al 2010, and Griffin et al 2013). The model then assigns phenotypic values to individuals, and fitness is derived. Figure two shows how the mean, skew, and kurtosis of the variance distributions changed for those 1500 generations.
"While the mean trait values change under selection, so also will the
G matrix, its orientation tending to shift in the direction of
selection. ... Genetic drift may also play a role in changing the G
matrix, but in this case the change will be random though on average
producing a proportional change in the constituent variances and