Substantive issues

Before getting into details of the choice of model type and coding, you need to make some decisions about what you want to model. As I understand the situation, you have epigenetic_age_acceleration (EAA)values for about 1700 individuals across 4 cities. Each individual has values for Smoking_Status, Sex and Age, and measurements of each of 5 heavy metal concentrations taken from near their residences. You want to model EAA as a function of those predictors.

First, both your ordinary least squares (OLS) model and your (incorrectly formulated) generalized additive model (GAM) suggest that you want to include interactions of each of the other predictors with the categorical predictor cities. It's not particularly clear what you hope to accomplish by that. You will not get a generalizable answer about the overall association of the other predictors with EAA. Rather, you will effectively get a set of 4 separate models, one for each value of cities. That might describe your data adequately but it's not clear how you would even think about extending those results to new situations.

To this non-expert, it seems that your model would be more useful if you used cities as a categorical predictor to allow for systematic overall differences in EAA among cities, but not to include interactions of cities with the other predictors unless you have strong theoretical reasons to expect that the value of cities per se will affect the associations of other predictors with EAA.

Second, both of your proposed models do not include any interactions among the predictors other than cities. They thus do not allow for combinations of heavy metals to have different associations with EAA than the sums of their individual associations, or for the effects of the heavy metals to depend on Age or Sex or Smoking_Status. That seems to be very limiting.

Before you continue, I recommend careful study of the principles of regression modeling, for example in Frank Harrell's Regression Modeling Strategies (RMS). Chapter 4 seems directly relevant to decisions about the complexity of the model that you might be able to fit with 1700 observations, and Chapter 2 discusses strategies for flexible modeling of your continuous predictors (heavy metal concentrations and Age). Some choices need to be made based on your understanding of the subject matter.

OLS or GAM

Once you have decided what you want to model, both OLS and GAM can allow for flexibly fitting nonlinear functions of continuous predictors. Chapter 7 of An Introduction to Statistical Learning (ISLR) outlines the approaches. Regression splines most naturally fit into an unpenalized OLS model. They can be thought of as a type of GAM, but a GAM is typically thought of as a more general approach allowing for several types of flexible modeling, with smoothness of the fits determined by penalization. Although ISLR emphasizes separate modeling of each continuous predictor in a GAM, the gam() function in the R mgcv package can accept smoothing terms that effectively incorporate interactions among predictors.

OLS with regression splines can work quite well. Frank Harrell notes:

As an aside, I have found little advantage of GAM's over parametric additive regression spline models, which give simpler formal tests and confidence intervals, and provide formulas for prediction.

OLS

Consider an extension of the following type of unpenalized OLS model. I recommend the rms package for this type of work, as its rcs() spline implementation has better default knot placement in my opinion than splines::ns() and (once you learn how to use it) the package provides a coherent framework for analysis, presentation and prediction.

library(rms)
olsModel <- ols(EAA ~ cities + Sex + rcs(Age) + Smoking_Status + rcs(Zn) + rcs(Cu) + rcs(Hg) + rcs(Mo) + rcs(Pb)

The default rcs() terms use 4 knots (3 coefficients to estimate) for each continuous predictor. If Smoking_Status is binary, then this model only requires estimating (beyond the intercept) 3 coefficients for cities, 1 for Sex, 3 for Age, 1 for Smoking_Status and 15 for all 5 heavy metals (3 each). In practice, you might want to consider working with log transformations of the heavy metal concentrations instead of their raw values, but still using splines for fitting them.

The above is the type of "saturated model" that Harrell recommends as a starting point in Section 4.1.1 of RMS.

With 1700 observations you could consider a much richer OLS model with interactions; you might envision a model with up to about 100 coefficients without overfitting. For example, including an interaction of Sex with each of the spline-fit heavy metals would only add 15 more predictors; the same for Smoking_Status if that's binary. In general, I find it safest to use "*" for interactions instead of trying to specify particular interactions with ":" as you did. That ensures that all lower-level terms are also included appropriately. The anova() function for rms objects (like those produced by ols()) provides a very helpful summary of overall significance of each predictor, including its nonlinear terms and interactions.

Sometimes interactions aren't needed beyond the linear components of the splines (specified by "%ia%" instead of "*" in rms syntax), substantially cutting down on the number of coefficients. Restricting to linear interactions could allow you to evaluate interactions among the heavy metals. Or, if you really want to do so despite my misgivings, to include interactions of cities with the heavy metals.

GAM

The main potential advantage that I see with a GAM of your data instead of OLS is if you want to allow for nonlinear interactions among the heavy metals. All of the simpler interactions of heavy metals with the other predictors, or even linear interactions among them, can be handled by OLS as noted above.

If that's what you want, heed this warning from Wayne:

I'd emphasize that GAMs are much more flexible than GLMs, and hence need more care in their use. With greater power comes greater responsibility... You don't have to understand exactly how GAMs work to use them, but you really need to think about your data, the problem at hand, your software's automated selection of parameters like smoother orders, your choices (what smoothers you specify, interactions, if a smoother is justified, etc), and the plausibility of your results.

If you want to use GAMs you have to devote some serious effort. It's not just a matter of how to code them; it's a matter of understanding enough about them to make intelligent choices.

Happily, there are resources to help with that. The help page for gam.models in mgcv deserves careful study. As best as I can tell, it seems to contain answers to all your questions.

For example, a model like this:

y ~ x1 + x2 + s(x3) + s(x4,x5)

treats x1 and x2 as unpenalized, uses a single s() smooth for x3, and a two-dimensional s() smooth combining x4 and x5, using "the default basis for the smooths (a thin plate regression spline basis for each), with automatic selection of the effective degrees of freedom for both smooths." That said, "The above assumes that x4 and x5 are naturally on similar scales... If this assumption is false then tensor product smoothing might be better (see te)."

If you want a "main effect plus interaction" structure:

Such models should be set up using ti terms in the model formula. For example:

y ~ ti(x) + ti(z) + ti(x,z)

For interactions of continuous and categorical predictors:

by variables are the means for... letting smooths ‘interact’ with factors or parametric terms... Note that when using factor by variables, centering constraints are applied to the smooths, which usually means that the by variable should be included as a parametric term, as well.

That's where your GAM interaction model was in error: you did not include cities as a parametric term (outside the smooth). I'm also not sure how including Hg as both a smoothed and as a parametric predictor would work.

If you need to go beyond that, there are nearly 1000 questions on this site tagged generalized-additive-model. The over 100 frequently viewed pages probably contain answers to other questions that you might have. In particular, Gavin Simpson is an expert in GAM who has contributed many helpful answers.