Calculating interactions for continuous variables with emmeans::emtrends()

[simRock]

I am trying to estimate the joint interaction for continuous variables with the emmeans::emtrends() function but I am having trouble doing so. Any help would be greatly appreciated it. See example below

library("tibble")
n=1e4
simdata <- tibble(
  age = rnorm(n, mean=30, sd=2),
  bmi = rnorm(n, mean=22, sd=2),
  age_c = age - mean(age),
  bmi_c = bmi - mean(bmi),
  y = rnorm(n, mean=2 + 2*age_c + 3*bmi_c + 4*age_c*bmi_c, sd=2)
)
model_y <- glm(y~age_c*bmi_c, family = gaussian(link = "identity"), data=simdata)
summary(model_y)
#> 
#> Call:
#> glm(formula = y ~ age_c * bmi_c, family = gaussian(link = "identity"), 
#>     data = simdata)
#> 
#> Deviance Residuals: 
#>    Min      1Q  Median      3Q     Max  
#> -7.773  -1.382   0.014   1.365   7.950  
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 1.990902   0.020204   98.54   <2e-16 ***
#> age_c       1.999852   0.010100  198.00   <2e-16 ***
#> bmi_c       3.000901   0.009994  300.27   <2e-16 ***
#> age_c:bmi_c 4.003130   0.005000  800.55   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 4.08059)
#> 
#>     Null deviance: 3217872  on 9999  degrees of freedom
#> Residual deviance:   40790  on 9996  degrees of freedom
#> AIC: 42447
#> 
#> Number of Fisher Scoring iterations: 2


library("emmeans")
#mean difference in y for 1 unit increase in age_c across levels of bmi_c
(beta_age_c_bmi_c0 = 2) 
#> [1] 2
(beta_age_c_bmi_c1 = 2 + 3.99)
#> [1] 5.99
emmeans::emtrends(model_y, 
         specs = "bmi_c",
         var = "age_c", 
         at = list(bmi_c= c(0, 1)))
#>  bmi_c age_c.trend      SE   df lower.CL upper.CL
#>      0       2.000 0.01010 9996    1.980    2.020
#>      1       6.003 0.01112 9996    5.981    6.025
#> 
#> Confidence level used: 0.95

#mean difference in y for 1 unit increase in bmi_c across levels of age_c
(beta_bmi_c_age_c0 = 3)
#> [1] 3
(beta_bmi_c_age_c1 = 3 + 3.99)
#> [1] 6.99
emmeans::emtrends(model_y, 
         specs = "age_c",
         var = "bmi_c", 
         at = list(age_c= c(0, 1)))
#>  age_c bmi_c.trend       SE   df lower.CL upper.CL
#>      0       3.001 0.009994 9996    2.981    3.020
#>      1       7.004 0.011227 9996    6.982    7.026
#> 
#> Confidence level used: 0.95

I want to get the mean difference in y for 1-unit increase in age_c AND 1-unit increase in bmi_c. I would like to obtain it using emtrends() but I am not sure what the command is. I am trying to get the output below

(beta_bmi_c1_age_c1 = 3 + 2 + 3.99)
#> [1] 8.99