Thank you both @Jose M. Martin-Garcia and @Henrik Schachner, your replies have been very helpful!

I am trying now to construct a matrix S, where the element ij contains the total population inside the previous circle, excluding the ones of cities i and j.

I have found that the populations can be obtained by using:

population = CityData[#, "Population"] & /@ cities

But I am not sure how to construct that matrix. Any idea would be much appreciated.

Hello David,

I'd suggest to use entities for the cities:

cities = {
   Entity["City", {"Istanbul", "Istanbul", "Turkey"}], 
   Entity["City", {"Moscow", "Moscow", "Russia"}], 
   Entity["City", {"London", "GreaterLondon", "UnitedKingdom"}], 
   Entity["City", {"Berlin", "Berlin", "Germany"}], 
   Entity["City", {"Madrid", "Madrid", "Spain"}]
}

Define also the two cities of reference:

madrid = Entity["City", {"Madrid", "Madrid", "Spain"}]
berlin = Entity["City", {"Berlin", "Berlin", "Germany"}]

Compute the distance between them (units depend on whether you are, so in the US we get miles):

In[]:= dist = GeoDistance[madrid, berlin]
Out[]= Quantity[1162.07, "Miles"]

Now you can query:

In[]:= Select[cities, GeoDistance[madrid, #] < dist &]
Out[]= {
   Entity["City", {"London", "GreaterLondon", "UnitedKingdom"}], 
   Entity["City", {"Madrid", "Madrid", "Spain"}]
}

This performs an independent call to WA per city, which will be slow for long lists. We can do better in those cases by fetching all positions at once in advance.