diff --git a/.github/workflows/ubuntu.yml b/.github/workflows/ubuntu.yml
index acd8ec0e3..352e3a145 100644
--- a/.github/workflows/ubuntu.yml
+++ b/.github/workflows/ubuntu.yml
@@ -133,7 +133,11 @@ jobs:
         run: |
           CPPFLAGS=$CPPFLAGS CC=$CC CXX=$CXX PYTHONPATH=. cmake examples/plugin
           make
-          PYTHONPATH=. python examples/plugin/test_plugin.py
+          PYTHONPATH=. python -We examples/plugin/test_plugin.py
+
+      - name: Run Python divergent optima example
+        run: |
+          PYTHONPATH=. python -We examples/gss_divergent_optima/gss_divergent_optima.py 100 0
 
       - name: Validate the sdist
         run: |
diff --git a/examples/gss_divergent_optima/gss_divergent_optima.py b/examples/gss_divergent_optima/gss_divergent_optima.py
index 3121ddc62..9aed27b0e 100644
--- a/examples/gss_divergent_optima/gss_divergent_optima.py
+++ b/examples/gss_divergent_optima/gss_divergent_optima.py
@@ -8,24 +8,42 @@
 N = int(sys.argv[1])
 rho = float(sys.argv[2])
 
-moving_optimum_deme_0 = fwdpy11.GSSmo(
-    [
+yaml = f"""
+time_units: generations
+demes:
+  - name: alpha
+    epochs:
+      - start_size: {N}
+  - name: beta
+    epochs:
+      - start_size: {N}
+migrations:
+  - demes: [alpha, beta]
+    rate: 1e-2
+"""
+
+demography = fwdpy11.ForwardDemesGraph.from_demes(
+    yaml, burnin=10 * N + 200, burnin_is_exact=True
+)
+
+moving_optimum_deme_0 = fwdpy11.GaussianStabilizingSelection.single_trait(
+    optima=[
         fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
         fwdpy11.Optimum(when=10 * N, optimum=1.0, VS=1.0),
-    ]
+    ],
 )
-moving_optimum_deme_1 = fwdpy11.GSSmo(
-    [
+moving_optimum_deme_1 = fwdpy11.GaussianStabilizingSelection.single_trait(
+    optima=[
         fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
         fwdpy11.Optimum(when=10 * N, optimum=-1.0, VS=1.0),
-    ]
+    ],
 )
 
 # The covariance matrix for effect sizes.
 # The marginal Gaussians will have mean zero and sd = 0.1
 # The deviates will be highly positively correlated.
 sd = 0.1
-covariance_matrix = np.matrix([0.999] * 4).reshape((2, 2))
+covariance_matrix = np.array([0.999] * 4).reshape((2, 2))
 np.fill_diagonal(covariance_matrix, 1)
 covariance_matrix *= np.power(sd, 2)
 
@@ -43,9 +61,7 @@
     ],
     "recregions": [fwdpy11.PoissonInterval(0, 1, rho / (4 * N))],
     "rates": (0, 1e-3, None),
-    "demography": fwdpy11.DiscreteDemography(
-        migmatrix=np.array([0.99, 0.01, 0.01, 0.99]).reshape((2, 2))
-    ),
+    "demography": demography,
     "simlen": 10 * N + 200,  # 200 gens past optimum shift
     # Specify one gvalue object per deme
     "gvalue": [
@@ -63,6 +79,7 @@
 pop = fwdpy11.DiploidPopulation([N, N], 1.0)
 rng = fwdpy11.GSLrng(1010)
 fwdpy11.evolvets(rng, pop, params, 100)
+assert pop.generation == 10 * N + 200
 md = np.array(pop.diploid_metadata, copy=False)
 df = pd.DataFrame.from_records(md[["deme", "g", "w"]])
 g = df.groupby(["deme"])
diff --git a/examples/plugin/test_plugin.py b/examples/plugin/test_plugin.py
index 1b9578dc3..75078469e 100644
--- a/examples/plugin/test_plugin.py
+++ b/examples/plugin/test_plugin.py
@@ -22,7 +22,12 @@
     "sregions": [fwdpy11.GaussianS(0, 1, 1, 0.25)],
     "recregions": [fwdpy11.PoissonInterval(0, 1, 5e-3)],
     "rates": (0.0, 5e-3, None),
-    "gvalue": fwdpy11.Additive(2.0, fwdpy11.GSS(VS=1.0, optimum=1.0)),
+    "gvalue": fwdpy11.Additive(
+        2.0,
+        fwdpy11.GaussianStabilizingSelection.single_trait(
+            optima=[fwdpy11.Optimum(VS=1.0, optimum=1.0, when=0)]
+        ),
+    ),
     "prune_selected": False,
 }